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A study has been made of the solution of the three-dimensional flow 
field for a flow-through nacelle. Both inviscid and viscous-inviscid 
interacting solutions were examined. Inviscid solutions were obtained 
with two different computational procedures for solving the three- 
dimensional Euler equations. The first procedure employs an allernating- 
di recti on -implicit numerical algorithm, and required the development of a 
complete computational model for the nacelle problem. The second 
computational technique employs a fourth-order Runge-Kutta numerical 
algorithm which was modified to fit the nacelle problem. Viscous effects 
on the flow field were evaluated *vi th a viscous-inviscid interacting 
computational model. This model was. constructed by coupling the explicit 
Euler solution procedure with a "lag-entrainment" boundary layer solution 
procedure in a global Iteration scheme. The computational techniques have 
been used to compute the flow field for a long-duct turbofan engine 
nacelle at free-stream Mach numbers of 0.80 and 0.94 and angles-of-attack 
of 0° and 4°. 

The numerical experiments show that_for predicting the flow inside 
the nacelle's duct, the viscous effects are extremely Important and both 
the external and internal boundary layers and wakes must be simulated. 

The internal boundary layer creates a pressure gradient in the nacelle's 

duct. The external boundary layer and its associated wake displace the 

streamlines away from the external surface. The displacement of the 
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streamlines effectively decambers the nacelle's airfoil and weakens the 
compression at the trailing. edge. This gives a less positive exit 
pressure and hence a less positive overall internal pressure level which 
agrees better with experiment than the inviscid computations. Therefore, 
if simulating the correct mass flow through a flow-through nacelle's duct 
is important, then viscous effects must be included in the computational 
model . 

In contrast to the internal surface, viscous effects were relatively 
unimportant for predicting the flow on the external surface of the 
nacelle. Good agreement is shown between the computational results of 
both Euler numerical procedures and wind-tunnel data on the external 
surface of the nacelle. The solutions exhibit the proper three- 
dimensional behavior at both angles of attack and correctly reflect the 
qualitative and quantitative results at both Mach numbers. 

The solutions of both Euler computational techniques exhibited a 
total pressure loss on the internal surface of the nacelle. An 
investigation of the loss proved that it was the result, of the flow 
stagnating on the external surface and expanding around the sharp 
discontinuity, in the surface of the nacelle at its leading edge. The 
studies indicate that -the use of C-type grids could probably eliminate the 
loss. However, for sharp leading edges or where an H-type grid is 
desirable. It appears that some (problem-dependent) total pressure loss is 
inherent in numerical Euler-equation solutions. 

Even though reasonably accurate engineering solutions were obtained 
with the Implicit computational procedure, a weak Instability was 
discovered in It when applied to the three-dimensional nacelle problem. 
This instability prevented the implicit solutions from actually converging 
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CHAPTER I 
INTRODUCTION 

Since fuel costs are projected to rise significantly over the next 
several decades, further growth of the air transportation system becomes 
increasingly dependent .on advances in aircraft technology and design. 1 
One important area of transport technology in which there is a potential 
for significant improvement is the integration of the propulsion system 
with the airframe. To tap this potential, engineers must increase their 
understanding of the aerodynamic interactions between the various- compo- 
nents of the propulsion system and airframe well beyond the understanding 
which now exists. 

Both experimental and theoretical research is required. 

Experimental studies of these interactions necessitate expensive, complex 
models to simulate adequately the inlet and exhaust flows. Consequently, 
experimental investigations are only practical for limited parametric 
studies. For analysis and design optimization, increasing attention is 
being given to the use of computational methods which have been validated 
by a few discrete experiments. For example, both the Airbus 310 and the 
Boeing 757 wings have been designed with the aid of numerical analysis. 2 

Because of the geometrical complexity of transport configurations 
and computer storage limitations, current computational design studies of 
this type have mainly been limited to inviscid potential -flow methods. 
These methods, which range from panel techniques to solutions of the full 
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transonic potential equation , 3 " 5 have proven to be useful tools in the 
aircraft design process# Nevertheless, they do not adequately simulate 
flows where rotational ity is important. 

To account for rotational ity, even in the Inviscid case, the 
application of the Euler equations is required.- Solving the inviscid, 
variable entropy, compressible flow equations allows the solution to 
capture strong shocks and simulate the jet exhaust flow without special 
treatments. In addition, vortex sheets may be captured, and as recently 
shown by Rizzi and Erickson , 6 a Kutta condition may not need to be 
explicitly enforced. The Euler equations also are thought to yield unique 
solutions, whereas the full potential equation can yield multiple 
solutions as shown by Salas^ and Steinhoff and Jameson . 6 Despite their 
greater potential, most solutions of the Euler equations have been either 
two-dimensional, or for relatively simple three-dimensional 
conf i-gurati ons , or on coarse grids .^"^ 3 Finally, strong interactions can 
occur between the viscous boundary layer and the main stream even when the 
boundary layer does not separate. Hence, the influence of viscosity on 
the flow field must also be accounted for in order to simulate the 
physically realistic case. 

Therefore, a study of the solution of the three-dimensional flow 
field for a flow-through nacelle using the Euler equations and a viscous- 
inviscld interacting computational model was initiated. A flow-through 
nacelle was selected for the study since.. the engine nacelle is a 
fundamental component of a transport's propulsion system. The objective 
of .this research was to investigate the flow field about the flow-through 
nacelle and the importance of simulating the viscous effects In a 
computational model for solving this problem. It included investigating 
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the suitability of the Euler equations as the basis .for the computational 
model and the problems associated with obtaining numerical solutions for 
th1s._type of configuration. In addition, it involved consideration of the 
advantages and disadvantages of the basic types of numerical algorithms 
and solution techniques suitable for flows of this complexity. In order 
to focus specifically on the application of computational solutions to the 
flow-through nacelle problem, only state-of-the-art algorithms were 
considered.. 

To conduct the study, two separate. Euler computational procedures 
were investigated, each of which appeared to offer certain distinct 
advantages. The first procedure employed an alternating-direction- 
implicit numerical algorithm and required development of a computational 
model specifically geared to the nacelle problem. An interim report on 
this phase of the research is given in reference 14. The second procedure 
investigated involved modification of an existing fourth-order Runge-Kutta 
algorithm to fit the nacelle problem. Viscous effects on the flow field 
were evaluated with a viscous-inviscid interacting computational model. 
This model was constructed by coupling the explicit Euler solution 
procedure with a "lag-entrainment" boundary layer solution procedure in a 
global iteration scheme. 

The computational techniques have been used to compute the flow 
field for a long-duct turbofan engine nacelle at free-stream Mach numbers 
of 0.80 and 0.94 and angles-of -attack of 0° and 4°. The results are 
compared with experimental data. New insight into the mechanics of the 
Interactions between the Internal and external flows,, gained during these 
investigations, is discussed. Problem areas, both general and algorithm 
dependent, are Identified and investigated. The numerical performance of 



CHAPTER II 


MATHEMATICAL DESCRIPTION OF THE PROBLEM 
Governing Flow Equations 

The Euler equations mathematically describe the physical laws 
governing the motion of an inviscid compressible fluid with variable 
entropy. In the present solution procedure, the three-dimensional time- 
dependent Euler equations are normalized and written in strong 
conservation form for a Cartesian coordinate system. If body forces are 
neglected, these time-dependent equations for mass, linear momentum, and 
energy can be expressed in vector notation as 



3f 

3x 



3h 

dz 


= 0 


where 


(1) 


p 


pu 

pu 


pu 2 + p 

pv 

I s 

puy 

pw 


puw 

pE 


u(pE + p) 


( 2 ) 
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pv 


pW 

puv 


pUW 

pv 2 + p 

■ . Jl = ■ 

pvw 

pvw 


pw 2 + p 

v(pE + p) 


w(pE + p) 


In these equations, u, v, and w are velocities in the physical 
coordinate system (coordinates x,y,z), p is the density, and p is the 
pressure. The total energy, E, is given by 

E * c v T + j (u 2 + v 2 + w 2 ) (3) 

where T is the temperature, and c y is the specific heat at constant 
volume. The equation of state 

p = pRT (4) 

where R is the gas constant, completes the system of_equations. 

The major interest in using computational fluid dynamics for 
propulsion integration studies is in predicting the steady, state flow for 
any given configuration. Therefore, the solutions to only steady flows 
are considered in this analysis. In this case, the total enthalpy, H, 
where 

H = E + £ (5) 

P 

does not vary throughout the flowfield of the flow-through nacelle; and 
the energy equation, the fifth equation of the set (1), could be replaced 
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by the condition of constant enthalpy. However, In general for propulsive 
flows the enthalpy is not constant due to the jet exhaust. Therefore, to 
be consistent with the ultimate objective of this study, the full energy 
equation was solved along with the continuity and momenta equations. 

Computational Domain and Grid System 

A sketch of the three-dimensional computational domain illustrating 
the nacelle and domain boundaries in both the physical and computational 
spaces is presented in figure 1. Three dimensionality is produced by 
rotating the vertical cross-section about the axis of the nacelle, thus 
generating a cylindrical domain (see part (a) of the figure). To minimize 
computer run time and storage requirements, symmetry is_assumed about the 
vertical plane, and only one-half of the cylindrical domain is computed. 
When transformed to the computational space, the coordinate system becomes 
a rectangular three-dimensional domain (see part (b)). The axis, which is 
singular, transforms into .a plane forming one side of the domain, and the 
nacelle surfaces transform into a common internal plane as illustrated in 
figure 1(b). 

In the computational domain, the grid system constructed for the 
dlscretised solution procedures is body fitted (grid lines coincide with 
the nacelle surface and other boundaries) in order to facilitate 
implementation of the boundary conditions. It is a sheared, H-type 
computational grid. Figure 2 presents a vertical cross section of the 
grid in the physical space, again Illustrating the nacelle geometry and 
the various boundaries. The grid mesh In the circumferential direction Is 
generated by rotating the vertical cross-section about the axis of the 
nacelle. The grid spacing in the physical space is geometrically 







OUTFLOW 

BOUNDARY 
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stretched away from the nacelle with grid. points being clustered near the 
nacelle surface and near the leading and trailing edges. The basic 
computations were made with 58 axial grid planes (30 along the nacelle), 

29 grid planes in. the radial direction,, and 11 in the circumferential 
direction. Additional, calculations were made with grids, refined in the 
axial direction:, one containing 68 axial grid planes (40 along the 
nacelle), and another with 115 axial grid planes (59 along the nacelle). 

Analytical Boundary Conditions 

For the mathematical description of a physical problem to be well 
posed, whether the partial differential equations are to be solved in 
closed form or numerically, the correct boundary conditions must be 
included. Further, in numerically solving a set of equations in any 
finite computational domain, boundary conditions arise from two different 
sources. First, analytical boundary conditions are necessary for the 
problem to be well posed, and second, numerical boundary conditions arise 
from a need to complete the differencing equations in the numerical algo- 
rithm. The analytical boundary conditions will be discussed in this 
chapter, while the numerical boundary conditions will be discussed in the 
chapters which present the computational procedures. 

The flow-through nacelle problem includes inflow, outflow, symmetry, 
far-field, and impermeable surface boundaries. For a well posed problem, 
only part of the flow variables can be specified at each boundary, with 
the number and combination of the variables depending on the type of 
boundary. Characteristic theory provides a guide to both the number and 
the form of the boundary conditions. Cline 15 describes the application of 
the theory for this purpose. The theory gives the same results for an 


invlscld fluid as the energy method of Oliger and Sundstrom,^ and was 
used as a guide In determining the boundary conditions for this work* 

Since the mainstream is essentially perpendicular to the Inflow and 
outflow boundaries, reference plane characteristics will be used to 
determine the number and form of these boundary conditions. Reference 
plane characteristics are an approximate form of the characteristic 
equations; however, they are exact if the cross derivative terms are 
zero. For the three-dimensional Euler equations, the constant-y-constant 
z reference plane characteristics can be written 


characteristic direction 


compatability relation 


p Ht “ *3 


dw _ , 
p St " *4 


dp 2 du _ 2 . . , . , . 

dt 3 3t - a + a* 2 + * 


it" “ pa BF “ " a ^? + ^ 



and a is the local speed of sound. For a flow with the mainstream in 
the x-di recti on, the t|> terms on the right hand side of these equations 
are very small and are treated as constants, and the equations become 
exact when the terms are identically zero. The characteristic direc- 
tions of equations (6), (7), and (8) represent the projection of the flow 
pathline and Mach cones on the x-t plane. Subject to the approximate 
nature of the reference plane characteristics, they comprise a set of five 
equations in five unknowns and completely describe the flow in the 
computational domain. 

Figure 3 presents a sketch depicting the projection of the 
characteristic curves on the x-t plane. For a boundary with Inflow In 
the positive x-di recti on, the shaded area represents the computational 
domain. Since three compatibility equations are associated with the 
curve dx/dt * u, equation (6), there are four characteristics entering 
the computational domain and one leaving it. The four characteristics 
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entering the domain carry no information whatsoever about the internal 
solution, and therefore they must be replaced by four boundary 
conditions. .The form of the. boundary conditions is obtained from, the 
characteristic leaving the computational domain, equation (8). .Since it 
originates from within the domain, both p and u cannot be prescribed 
because they must obey the compatabi.l ity relation holding along the 
characteristic curve. Usually the total pressure, total temperature, and 
flow angle are supplied at an inflow boundary, or, an alternate 
combination consisting of a function of the total entropy, the static 
temperature, and velocities. 

For a boundary with outflow in the positive x-direction, the 
unshaded area of figure 3 represents the computional domain. The figure 
illustrates that there is only one characteristic entering the computa- 
tional domain, indicating one boundary condition must be supplied. It Is 
usually the static pressure, since, when combined with the previous inflow 
boundary conditions, it completely defines the free stream. 

The far field boundary can be either an inflow or an outflow 
boundary. Therefore the appropiate conditions described above are 
applied, 

The expression for the boundary conditions at the symmetry plane is 

if ■ 0 do 

where n_ is the normal to the symmetry plane. 

Finally, the surface boundary condition Is that the flow Is tangent 
to the wall, and is expressed by the equation 


-mpyir'ir • - - 
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n»V. = 0 (H) 

where V is the velocity at the surface and ^represents the unit vector 
normal to the surface. 

As mentioned above, in addition to these analytical boundary 
conditions, numerical boundary conditions necessary to complete the 
difference equations are required when actually carrying out the numerical 
process. They will be discussed in the next two chapters. 

Nacelle configuration 

A flow-through nacelle was selected for the study since the engine 
nacelle is a fundamental component of a transport aircraft's propulsion 
system. The particular flow-through nacelle chosen is depicted in figures 
1 and 2. The external surface of the nacelle consisted of a NACA 1-70-100 
inlet, a cylindrical section, and a circular arc afterbody. The Internal 
duct of the nacelle was cylindrical. For the computations, the leadin.g 
edge was made sharp as a simplification to the geometry. This particular 
configuration was chosen because it resembled, long-duct turbofan engine 
nacelles being proposed for current jet transports, and because experimen- 
tal-data for the isolated nacelle were available for comparison with the 
computational results. 


CHAPTER III 


IMPLICIT COMPUTATIONAL PROCEDURE 

The implicit procedure used for the flow-through nacelle problem 
employs the approximately-factored, alternating-direction implicit 
algorithm of Beam and Warming . 17 This algorithm has been applied success- 
fully to a number of two- and three-dimensional problems . 18 * 19 The 
principal advantage offered by implicit methods is that, if properly 
formulated, they theoretically have no stability limitations on the size 
of the time step when integrating the set of flow equations. Thus,, for 
obtaining steady-state solutions, fewer integration steps may be needed. 
The use of a large time step can in some cases further act to accelerate 
the convergence rate in a manner similar to relaxation schemes for 
elliptic problems. The principal disadvantage of implicit methods is the 
requirement for solving large sets of simultaneous algebraic equations. 
Thus while fewer time steps may be required than for explicit methods, 
more computational work per time step is usually needed. However, the 
previous applications of the Beam and Warming. algorithm have indicated 
favorable Improvements in overall computational efficiency. The applica- 
tion of this algorithm to solving the three-dimensional Euler, equations 
for flow-through nacelles required .the development of a computational 
model specifically geared to the problem. This chapter describes that 
computational model. 
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Transformed Euler Equations 

Equation (1) presented the Euler equations written for a Cartesian 

coordinate system. In. .the present finite difference technique, the 

equations are transformed to the computational space before they are 
integrated. The transformed version of the three-dimensional Euler 
equations can be written 


A A 



( 12 ) 


where t is the normalized time, 5 , 5, and n are the spatial coordi 
nates in the transformed plane, and 


q = J 4 

F * J(c x l + ?y£ + 5 z h) 

6 - J(5 x f + + 5 z h) 

H = 0(n x f + n y< £ + n 2< h) 


(13) 


where tb: vectors, _cj_, f_, j,, and h_ are the vectors for. the Euler 
equations written in Cartesian coordinates (see Chapter II). In the 
relations above, J .is the Jacobian of the transformation .from the 
physical plane (coordinates x,y,z) to the computational plane 
(coordinates ;, c, n ). The expression for the Jacobian Is 


19 


Vc z n + x cVn + Vs z 5 " x n y c z ; " 


xv z . - x.y z 
r n C 5 s n 


(14) 


Equations (12), (13), and (14) are valid for transformation between 
any two coordinate systems. In this particular case, the coordinates, 
x, y, and z, represent respectively the axial, vertical, and 
horizontal Cartesian coordinates in which the Euler equations were first 
written. The coordinates, c, and n, represent respectively 
functions of the axial, radial, and circumferential directions of the 
cylindrical computational domain with the stretched grid spacing. They 
can be expressed as 


t - f l(x> 

C ■ f 2 < r > 

n = f 3 (6) 

where 

r * (y^ + z^)^ 
e = tan ' 1 (y) 


(15) 


(16) 


and fj, f 2 * and f 3 are stretching functions. The ?, 5, and n 
coordinates "unwrap", to form the computational space with equal spacing 
between coordinate lines. A point to point correlation between the 




physical space and the computational space is given in figure 1 and 


discussed in Chapter II. 


Numerical Method 


Algorithm . The computational technique employs the approximately- 
factored alternating-direction implicit numerical algorithm of Beam and 
Warming in its delta form. 17 The delta form of the Beam and. Warming 
scheme has the advantage that the steady state solution, if one exists, is 
independent of the size of the time step. The algorithm, including 
stability and accuracy limitations, is described in detail for two dimen- 
sions by Warming and Beam;^ therefore only a brief outline of the method 
will be presented here. 

The. algorithm is developed by first applying an. Euler implicit 

A n+1 

formula between time levels n and n+1 to express the vector q of 

A n 

the transformed Euler equations (12) in terms of q . Then .by applying a 

A A A 

linearization procedure to the vectors F, G, and H, using a local 
Taylor expansion about q n , the Euler equations can be written in the 


I 1 M!?*" + lr B " J k c ")K +1 ‘ 


f " + !? 5n ♦. fc H") + - °(“ 2 > 


where I is the identity matrix. The term in the braces on the left-hand 

A n+1 

side of the equation (17) operates on Aq where 
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*n+l -n+1 «n 
Aq * q - q 


( 18 ) 


The terms A n , B n , and C n are Jacobian matrices defined as 


A" = B n » 4^ C n = 


aq 


aq 


aq 


(19) 


Equation (15) is then approximately factored giving 


(' +it lr An )( I + it lr ,,n )( I + it i c >5 


n+1 


-“(It ^ It s " ♦ It 5 ") 


( 20 ) 


The factorization of equation (17) produces additional terms which are on 
2 

the order of At and hence do not destroy the temporal accuracy of the 

scheme. However the present studies indicate that, in three dimensions, 

the additional terms may seriously degrade, the stability characteristics 

of the algorithm. This property of the algorithm will be discussed in 

detail later. The spatial derivatives In equation (20) are approximated 

by second-order central -difference operators yielding a block-tri diagonal 

system of equations to be solved. in each of the three directions 

C, £, and n. The alternating-direction-implicit sequence for 

A n+1 

determining q is 



r were added to the appropriate factors on the Implicit left-hand side. The 

I" ' 

jp value of the explicit damping coefficient, u> e , Is set equal to at, 

?• ] and the value of the Implicit damping coefficient, , Is set equal 

If* to 2w e In the manner of Pulliam and Steger.^ For the numerical 

[>• calculations, the derivatives In the dissipation terms are approximated by 


_-% m t 

\ ' v> • • • \ 
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finite-difference formulas. It should be noted that although numerical 
dissipation can stabilize a solution process, it has the disadvantage of 
causing a loss of. accuracy in regions of strong gradients. 

Recently, Abarbanel , Dwoyer, and Gottlieb 22 have proven that the 
undamped Beam and Warming scheme in three dimensions is weakly, but uncon- 
ditionally, unstable for the Euler equations. Their work, which has been 
accomplished since the present research was initiated, shows that the 
instability is very weak and is only present in the very long wave- 
lengths. .Numerical experiments conducted in the present investigation 
indicate that the dissipation terms described above do not fully stabilise 
the solution process. However, the experiments also show that the weak 
long-wavelength instability takes a very large number of time steps to 
develop near the body. For a reasonable number of time steps (~ 1000), 
reasonably accurate solutions were obtained. 

Metric calculation . If the scheme is applied in a uniform free 
stream, it is expected.. that uniform free-stream conditions would be 
exactly maintained. In three dimensions,, however, errors in the uniform 
flow can be introduced since the transformed flux terms being differenced 
In equation (20) contain the metric derivatives. These errors can be 
avoided by using the proper averaging technique when numerically 
calculating the metric derivatives as pointed out by Pulliam and Steger. 2< 

Since three-point central -difference operators are used to evaluate 
the flux terms of equation (18), the errors will exactly cancel if, for 
example, c is calculated by the equation 
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U x>1,j,k * T 1tj>k lC(6 j y 1.k+l + yi,k-l )/2] - 

* ^ 6 k z i,j+l + 6 k z 1,j-l^ 2 ^ ~ ^ 6 k y 1,j+l 
+ 6 k y i,j-l )/2]C(6 j z i,k+l " 6 j z i,k-l )/2]} (24) 


and the other metric derivatives are calculated in a similar manner. In 
equation (24), &. and 6 k are the three point central difference 

operators in the j and k directions,, respectively. The metric 
derivatives were computed by this technique to hold numerical errors to a 
minimum and consequently enhance the stability of the computational 
procedure. 

Boundary Conditions 

Inflow boundary . The analytical boundary conditions were covered in 
the previous chapter. The numerical boundary conditions, necessary. to 
complete the differencing equations in the finite difference scheme, will 
now be specified for the implicit algorithm. At the inflow boundary, the 
analytical boundary conditions are met by specifying the total pressure, 
total temperature, and inflow angle. The static pressure is tn« remaining 
unknown on the boundary and must be numerically determined. It must obey 
the compatabillty relations along the characteristic connecting the Inflow 
boundary with the Interior solution. Comparability .is approximated by 
extrapolati ng t he pressure from within the computational domain. 

Outflow boundary . At the outflow boundary, the analytical boundary 
condition is. the specification of the static pressure. Two basic methods 
were used to specify the static pressure for the implicit numerical 


scheme. The first -technique evaluated was setting the boundary pressure 
equal to the free-stream value, l.e. 


In addition to this method, a radiation boundary condition based on a 
solution to the three-dimensional wave equation 



where 


( 26 ) 


was tested. Bayliss and Turkel 24 derived the radiation boundary condition 
and applied it in two dimensions with good results. During the present 
research, the radiation boundary condition was transformed to a three- 
dimensional form and applied to the flow-through nacelle problem. 

In the far field where pertubations are small, the three-dimensional 
Euler equations assume the form of equation (26) when linearized and 
transformed by the relations 



If spherical coordinates are also introduced, 


* 2-2 2 2 

r* * c + y + z 

6 = tan '■ 

$ = tan" 1 ^-|-j 

equation (26) has solutions of the asymptotic form 


(29) 


p = til - [y». 1) 
r 


(30) 


A 

for large x and r. The functional form in equation (30) represents an 
outgoing spherical wave solution to (26). The specific function, f, is 
not usually known. However, the condition 


8 3r r 


(31) 


is exact for ail functions which identically have the form (30). 
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Equation .(31) Is. the analytic boundary condition used at the outflow 
boundary. Transforming back to the physical plane (coordinates t,x,y,z) 
gives 



and the superscript, n, indicates the time step. During transient 
periods, equation (32) allows impinging pressure waves to pass through the 
outflow boundary more effectively than equation (25). At steady state, 
equation (32) reduces to equation (25). In the implementation of both 
boundary equations, (32) and (25), the flow quantities other than the 
static pressure are obtained from the interior of the computational domain 
by zeroth order extrapolation. 

In addition to the two basic outflow boundary conditions just 
described, the nonreflecting outflow boundary condition of Rudy^ was 
tested briefly with_m1xed results. The Riemann-lnvariant method described 
in the next chapter was also tried in conjunction with the implicit 
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computational technique. An equally brief Investigation of this 
combination failed to yield a satisfactory solution. 

Far-Field boundary . The treatment of the far-field boundary depended 
on whether outflow boundary condition (25) or (32) was used. When outflow 
boundary equation (25) was used, all the flow quantities on the boundary 
were obtained by zeroth order extrapolation from interior points. When 
the radiation equation, (32), was applied, it was simultaneously applied 
at both the outflow and far-field boundaries (see figure 2) as suggested 
by Bayliss and Turkell. 2 ^ In the application of the radiation condition 
at the far-field boundary, p^ in the right-hand side of equation (32) 
was replaced by the local value of the pressure at the previous time step. 

A more precise treatment of the far-field boundary would.be to test 
each point on the boundary at each time step for an inflow or an outflow 
condition and treat the point accordingly. However, since the boundary is. 
relatively far from the nacelle, and the velocities normal to the boundary 
are very low in comparison to the free-stream velocity, the less 
complicated treatments were chosen for the calculation at 0° angle of 
attack. For calculations with the nacelle at angle of attack, inflow 
boundary conditions are assigned to the lower half of the far-field 
boundary and outflow boundary conditions to the upper half. 


i 


i 


Surface boundary . For an In viscid fluid, the boundary condition at 
an impermeable surface Js that the flow is tangent to the surface. This 
condition is expressed at the nacelle surface by. equating the 
contravariant velocity which is not tangent to the surface to zero. 


I ■ 6 U + ij + iJL * 0 

A Jr 4 


( 34 ) 
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By using this expression .and combining the three momentum equations In the 
manner of Pulliam and Steger, 22 

E (x-momentum equation) 

A 

+ £(y -momentum equation) 

+ E z (z-momentum equation) = 0 (35) 

the following relation which is independent of time can be obtained for 
the surface pressure: 

ap ? + bp^ + cp^ - d (36) 

where 


a s Ec + £ c + Ec 
*x^x *y v y s z v z 

2 2 2 
b * < * «; + 

c “ ^x ? x + + ^z ? z 

d = *Pii(C x u ; + S y v ; + C z w ; ) 

-pM(5 x u n + y n + S 2 w n ) 

and y and Jrf are the contravariant velocity components tangent to the 




nacelle surface. 


U. = c x £ + CyV. + C 2 w 

(38) 

W = n x< u + n y v. + n z w 

Equation (36) is solved for the surface pressure .by first 

extrapolating p, U, and W_ along the £ coordinate to the nacelle 

using a first-order procedure. Then p^ is expressed as a second-order 

one-sided finite-difference formula, and p, and p are written as 

c n 

second-order central difference formulas. The resulting equation is then 
approximately factored, and a series of tridiagonal equations in 
the ? and n directions are numerically solved to yield the surface 
values. 

Leading and trailing edges . The implicit finite-difference computa- 
tional procedure requires a direct treatment of the leading and trailing 
edae boundaries, and a wide variety of the treatments was attempted in the 
present work. Solutions to the Euler equations at such points require 
careful treatment in order to overcome the mathematical difficulties while 
at the same time maintaining correct physical behavior. For example, the 
condition of tangency of the flow to the nacelle surface necessitates that 
stagnation conditions exist at the sharp leading and trailing edges. 
However, for the grid spaclngs investigated, specifying stagnation 
conditions at the leading and trailing edges resulted in large jumps in 
the flow quantities in the immediate vicinity of these points. A more 
accurate solution is obtained by the following approximate treatment of 
the leading and trailing edge boundaries. 
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At the leading edge, the density and contavariant velocities are 
extrapolated along the internal and external surfaces of the nacelle. The 
condition that the total enthalpy .is constant plus the extrapolated, 
quantities determines .the thermodynamic and kinetic. properties ..of the 
leading edge flow. The internal and external surfaces are treated 
se para tely and the tangency condition is maintained on each surface. 

The densities and contravariant velocities are also extrapolated 
along the internal. and external surfaces to the trailing edge. The 
pressure at the trailing edge is . determined from the external values and 
the condition that the total enthalpy is. constant. A discontinuity in the 
density and in the magnitude, of the velocities is allowed between the 

internal and external surfaces. 

Best agreement between the inviscid computations and the 
experimental pressures on the inside of the nacelle were obtained when a 
Kutta-like condition was adopted at the trailing edge. It consisted of 
setting the flow angle at the trailing edge equal to the angle of the 
internal surface of the nacelle. The primary results presented for the 
implicit procedure in Chapter VI are calculated using this boundary 
condition. However, further numerical studies show that instead of giving 
the correct Inviscid solution, the Kutta-like condition actually models 
viscous effects present in the experimental data. These results and their 
implications are presented in Chapter VII. 






CHAPTER IV 


EXPLICIT COMPUTATIONAL PROCEDURE 

In. addition to the implicit finite-difference numerical procedure 
described in Chapter III, an explicit finite-volume procedure which 
employed a fourth-order Runge-Kutta numerical algorithm was also evalu- 
ated. The algorithm has been applied to both two and three-dimensional 
problems, and appears to give accurate and stable solutions in both 
cases . 11 " 13 The main advantage of explicit schemes is that the updated 
solution at a new time step is independently calculated at each grid point 
or cell. Thus, sets of algebraic equations do not have to be solved as 
they do for implicit methods. The disadvantage of explicit numerical 
algorithms is that the grid spacing imposes stability limitations on the 
size of the time step when integrating the flow equations. Thus, less 
computational work is required at each time step, but more time steps are 
normally needed. The limitation on the size of. the time step Is most 
restrictive for viscous solutions where the grid spacing must be very fine 
in order to adequately define the boundary layer. For inviscid solutions, 
the restriction is much less severe. . 

Unlike the implicit computational procedure, the explicit procedure 
was not completely developed in the present study, but was adapted to the 
nacelle problem from an existing computer code written by Jameson and 
Baker . 13 The adaptation consisted mainly of changing the logic of the 
code to allow for both internal and external hacell e surfaces on an H-type 
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grid. Essentially, the explicit interior point algorithm and boundary 
treatments were inserted into the overall framework developed for the 
implicit code. The adaptation also involved reprogramming some of the 
boundary conditions for the H grid. This chapter describes the applica- 
tion of the explicit algorithm to the nacelle problem. 

Finite Volume Formulation . 

The basic principles of the explicit procedure are covered, 
thoroughly in references li, 12, and 13; therefore only a brief outline of 
the procedure will be presented here. The Euler equations, (1), can be 


written in integral 

form as 



3 i 

r r 




a~c' J 

1 qdV + / 

T • ds =0 


(39) 

J 

v J s 




where q 

is the vector of conserved quantities (dependent variables) 

presented in equations (1), V 

is the 

volume of the domain, S is its 

surface ■ 

area, and 1 

P is 




"pu. 

PV, 

p w 



pu 2 + p. 

puv, 

puw 


F * 

PUV, 

PV 2 + p. 

pvw 

(40) 


puw, 

pvw, 

pw 2 + p 

• 


.puH, 

pvH, 

pwH 


The commas separate 

the x, y, 

and z, 

, or physical , components of F. 


No transformations are necessary with the finite-volume formulation. 


r- r 


Numerical Method 


Algorithm . To integrate the equations over the computational domain, 
the explicit procedure uses the four-stage Runge-Kutta numerical algorithm 
in its standard form. For a linear system of equations 

ff*P(q)=0 (41) 

the four-stage scheme can be written as 

q(0) = q n 

qd) . q (0) . ai p(0) 

q(2) - q(0) _ a ^p(l) 

(42) 

q(3) = q (0) . ^(2) 
q(4) a q (0) . p(3) 
q n+ l = q(4) 

where q n and q n+ * are the values of the dependent variables at time 
steps n and n+1 respectively. For the standard scheme, the values for 
the coefficients, a, are 
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These coefficients give the scheme a numerical stability bound 
which allows a maximum Courant number of 2V?. 


Before the Euler equations can be numerically integrated by using 
the Runge-Kutta algorithm, they must be.discretised. In the present 
procedure, dissipation terms are also appended to the discretised 
equations before the solution process is begun. The next two sections 
describe the discretization process and dissipation terms. 

Discretised equations . The Euler equations are spatially discretised 
by approximating them in integral form in each cell of the computational 
domain. Note that the computational grid divides the domain into a system 
of adjoining hexagonal cells. The discretised equations for a cell are 

= ° < 44 > 

where the indices i,j,k identify the center of the cell, and v is the 
differential cell volume which is independent of time. The summation on 
the subscript, *, denotes summation on all six faces of the hexagonal, 
cell. At each face, the dot product between the face area, S, and the. 
flux tensor, F, is evaluated as 


S»F = S V F 


-t-Sy^y + s z F z 


( 45 ) 


where the subscripts x, y, and z Indicate components in the physical 
directions. Since the intermediate solution, the solution at time level 
n, is known at each cell Center, the flux at each cell face is evaluated 
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by averaging the fluxes at the centers of adjacent cells. For example, 
for. a typical face, I, located at 1+l/2,j,k with adjacent cell centers 
at i,j,k and 1+1, j,k, the mass flux at time level n would be 

( pu) I ■ 7 I (pu )i ,j ,k + < 46 > 

This method of averaging amounts to using central differences in a finite 
difference formulation, and yields second order spatial accuracy in the 
absence of any dissipation terms. 

Approximating the time derivatives with a forward difference formula 
and solving for q n+ * gives 


q 


n+1 




(47) 


which completes the discretization of the integral form of the Euler 
equations. 

Numerical dissipation . .Equation (47) is in a form compatible with 
the Runge-Kutta numerical algorithm. However, to insure numerical 
stability, dissipation terms are appended to the right-hand side of the 
equation before the solution process begins. The dissipation terms are 
the same as those used by Jameson et. al . in reference 13. The augmented 
form of the equations is 


q 


n+1 
i , j ,k 



At 

i * j >k 






(48) 




where represents thejilssipation terms In all three computational 

directions (directions 5, and n). 

D i,j,k “ (D ; + D c + D n ) 1 »j ,k ( 49 ) 

The damping terms are composed of third order terms (fourth-order 
differences) which prevent odd -even point decoupling , and first order 
terms (second-order differences) which prevent preshock oscillations. 

The third order terms are similar to the ones appended to the 
implicit procedure, and are formed by taking the fourth-order differences 
of the dependent variables along each of the computational coordinate 
directions. Taking the 5 direction as an example gives 

^1,j,k = 4+1/2, j,k * d i-l/2,j ,k ^ 5 °) 


where 


1 1+l/2,j,k 


. u (4) (v i : a 


*1+2 ’ 3q i+l + ^1 ’ q i-l ) j,k (5r 


and d^ 4 ^_j/2,j,k * s calculated by a corresponding formula. The 
term, (v^ + v^ + ^ ) j ^ ^ / 2At insures that the units are consistent in 
equation (48). Similar terms are formed for the 5 and n directions. 
These dissipation terms are spatially accurate to the third order, and 
hence do not compromise the accuracy of the spatial discretization of the 
Euler equations which is second order. 



In addition to the third-order dissipation just discussed, first 


order dissipation is necessary in the vicinity of shocks to prevent 
preshock oscillations. These terms are. formed by taking second-order 
differences of the dependent variables, and have the form 




e i +1/2 \t /2 ( q i+l " V 


i-1/2 _ x 

e i -1/2 At (q i " q i-l ) 


(52) 


where v and e form a switch sensitive to the second differences in 
pressure. 13 Their formulas for the face located at i +1/2 are 


v 


i.j.k 




(53) 


and 


=1+1/2 m ,ax < v 1+l,j,k* < 54 > 

Similar formulas are used to calculate v and e at 1-l/2,j,k. Near 
pressure discontinues such as. shocks, the first order dissipation terms 
become very large, and the scheme is first order accurate. Whenever the 
pressure variations are smooth, the first order terms are negligible, and 
the scheme remains spatially second order accurate. 
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The numerical dissipation just discussed is added to. equation (47) 
at all four stages of the Runge-Kutta algorithm. However, the damping 
terms are updated only at the beginning of each time step, and are not 
recalculated between the stages of the scheme. 

Implicit smoothing . In the present explicit method, like reference 
13, the stability bound is increased by applying implicit smoothing to the 
residuals which are calculated explicitly. The smoothing method Is 

pc 

similar to the class of numerical algorithms suggested by Larat, and 
takes the form 


(1 -5 «‘H1 


e " £ 6 n ) p i ,k = P i ,J,k 


(55) 


where p is the smoothed residual. The smoothing operators are applied 

A. 

as factors operating on P in the transformed directions. Therefore only 

% 

three separate tri diagonal equations must be solved to determine the 
smoothed value of . the residuals. 

Convergence acceleration . Convergence to a steady state solution is 

io 

accelerated in the computational procedure by two methods. First the 
maximum allowable time step is used for a given Courant number at each 
individual cell. Using a variable time step destroys the temporal 
accuracy of the solution technique, but that is of relatively little con- 
cern since the steady state Is the solution of Interest. A variable time 
step is equivalent to making the following modification to the Euler equa- 
tions 


21 ♦ ei 21 + ££ + IS. 

at ax ay az 


0 


*> '*5r 


(56) 
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where e Is a scalar multiplier and 1 is the Identity matrix. Thus 6 

can be chosen Individually for each cell to advance the solution at the 
maximum Courant number# 

The second method of enhancing convergence Is by including a forcing 
function in the Euler equations proportional to the difference between the 
local enthalpy and the free-stream enthalpy. With the addition of this 
term, the modified form of the equations becomes 


ia + it 

at ax 


+ !£.+ !£ L + aVu 

ay az 


H J 


( 57 ) 


where 


A(H - HJ 


"ap(H - H J 
apu(H - HJ 
apv(H - HJ 
apw(H - HJ 
.5 (H - H J 


( 58 ) 


and a Is a constant. At steady state. H • H. , so the steady state 
Euler equations remain unaltered when the enthalpy damping term is 
Included. Jameson also shows that In the absence of shocks, the flow 
would remain Irrotatlonal and homentroplc. 12 The enthalpy damping 
function acts similar to the term of the telegraph equation 


♦tt + * *xx + *yy 


<V 


5 *** ' 


( 59 ) 
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where In relaxation methods, plays an Important role In determ ining 
the rate of convergence. 

Boundrary Conditions 

Inflow boundary . The treatment of the inflow, outflow, and far-field 
boundaries is based on the Rlemann invariants for a one-dimensional flow 
normal to the boundary. ^ For a one-dimensional flow, the Riemann invari- 
ants and the characteristics along which they apply are 


characteri Stic Riemann invariant 



11 

C 

♦ -&T 

y - 1 

dx 


2a 

cit S u - a 

R = u 

"T^T 


( 60 ) 


where R is the value of the Riemann invariant, u is the velocity 
normal to the boundary, and a is the speed of sound. 

In applying the Riemann invariants to calculate the numerical 
boundary conditions, the values of u and a are determined from the 
flow variables in the regions where the respective characteristics origi- 
nate. For the inflow boundary depicted in figure 1, the characteristic 
originating from outside the computational domain is dx/dt * u + a. 
Therefore, the value of the corresponding Riemann Invariant Is calculated 
by the formula 


R 


U • h 



( 61 ) 
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where Is the unit normal to the -boundary and is the velocity of 
the f ree-st ream. - The Rlemann Invariant for the characteristic originating 
inside the computational region, dx/dt * u - a, is 

2a, 


int = — int* — ' T^T 


• n - 


‘int 


VS- *7 (r - + R 1nt> 


(62) 


where the subscript, .int, indicates the interior of the domain. In 
making this calculation, the values of Jij rit and a^ nt are computed from 
the properties of the flow variables in the cell next to the boundary. 

The equations for the two invariants, (61) and (62), can be. combined to 
yield the normal velocity and the speed of sound at the boundary. 


a b = 


Y - 1 

3 






(63) 


where the subscript, b, indicates boundary values. 

This procedure for calculating the normal velocity and speed of 
sound is equivalent to specifying one boundary condition and numerically 
calculating another. Consistency with the analytical boundary conditions 
discussed in Chapter III must be maintained. . Therefore, at the Inflow 
boundary where four-conditions must be supplied, the tangential velocities 
and entropy are set equal to their free stream values. The three velocity 
components, the speed of sound, and the entropy completely determine the 
dependent variables at the Inflow boundary. 
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Outflow boundary . The outflow boundary is treated in a similar 
manner. However, the characteristics entering and leaving the 
computational domain are reversed, resulting in the new equations for the 
normal velocity and speed of sound 

V* = I (R 1nt + R « ) 

(64) 

a b = Y 4 ' R 1nt ' R «' 

On the outflow boundary, the entropy and tangential velocities are 
extrapolated from the interior of the computational domain. Thus, the 
treatment of the outflow boundary. is also consistent with the number and 
form of the analytical boundary conditions given by characteristic theory. 

Far-field boundary . The far-field boundary is treated as an inflow 
boundary where the normal velocity is the vector sum of the velocities in 
the y and z directions, or v and w respectively. 

Surface boundary . The surface boundary treatment is very similar to 
the treatment for the implicit procedure. The concept is the same; 
however it_i.s applied slightly differently to the finite volume 
formulation used in the explicit numerical procedure. For the finite 
volume formulation, the pressure is the only flow quantity required on the 
body surface. A brief overview of the method, which was proposed_by 
Rizzi 27 and used in reference 13, follows. 

The boundary condition for a solid body in inviscld flow is that the 
flow is tangent to the surface, equation (11). Taking the total 
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derivative of the flow t 3 ngency condition with respect to time 


- + V_*grad) {V?n) ~ 0 


(65) 


(where V_ is the velocity at the surface and n_ represents the unit 
vector normal to the surface) and substituting it into the inner product 
of n and the momentum equations yields 


P V/(v»grad)n_ = n_*grad(p) 


( 66 ) 


for a stationary body. For a general nonorthogonal grid, equation (66) 
can be expressed in the form 




= q m2 ie_ 
9 


where 


(67) 


dXj 

*1 = dT 


ul _ ds 1 u 3 = d$ m2 „ 

W - W - » 9 ax. 


dt 


m2 _ i r iC 

3x„ 9x„ 


i i 


( 68 ) 


and double indices indicate summation as in tensor notation. 

In equation (67), Rizzi substitutes the. flow quantities at the 
nearest field point for those on the body surface. This yields a first 
order accurate boundary method .which is consistent with the accuracy of 


the second-order... interior scheme and gives an overall solution accuracy of 
second order. Thus, the gradient of the pressure along the coordinate 
line intersecting the surface, C, can be evaluated, and then the 
pressure can be extrapolated to the body. 


CHAPTER V 


VISCOUS-INVISCIO INTERACTING COMPUTATIONAL MODEL 

Chapters III and IV respectively describe an implicit and an 
explicit computational procedure for solving the three-dimensional Euler 
equations for the flow past a flow-through nacelle. However, the Euler 
equations, which model compressibility and rotationality , do not model 
viscous stresses. In the physically realistic case, strong interactions 
often occur between the viscous boundary. 1 ayer and the main stream even 
when the boundary layer does not separate. In these instances, modeling 
the frictional forces becomes essential if accuracy is to be maintained. 
Therefore, to complete the study of the flow field surrounding the flow- 
through nacelle, a procedure with which to assess the viscous effects was 
needed . 

To obtain a computational technique which simulated viscous effects, 
the explicit Euler solution procedure was coupled with a boundary layer 
solution procedure. The resulting viscous-inviscid interacting 
computational model is based .on a global iteration between the integration 
of the Euler equations and the boundary layer equations. The present 
chapter describes this interacting computational model. 

Boundary Layer Equations 

Since the objective of the viscous-inviscid interaction phase of the 
research was to evaluate the viscous effects on the nacelle pressures, a 


46 


47 


boundary layer solution procedure that was well validated was desired. 

The inviscid calculations and the wind-tunnel data indicate that the flow 
on the nacelle remains attached. Therefore an integral technique used in 
the direct mode was considered applicable.^® Green’s^ compressible, 
axi symmetric, "lag entrai nment" method solved in the direct mode was 
chosen because of its reliability. , 

The method involves the integration of three ordinary differential 
equations: the momentum integral equation, an entrainment equation, and a 

rate equation for the entrainment coefficient. 





1 + .075M 


2 (1 * .M.) \ 

(1 + .1M Z )l 


(69) 


The momentum integral equation,, the first equation, is obtaind by 
integrating in the direction normal to the wall both the continuity and 
streamwise momentum equations and combining the results. The entrainment 
equation, which is the second one, is obtained by integrating the 
continuity equation in the direction normal to tite_wall . The rate 
equation for the entrainment coefficient comes from a similar integration 
of the energy equation, and represents explicitly the balance between the 
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convection, production, diffusion, and dissipation of kinetic energy. 

In the previous equations, e is the momentum thickness and R the 
radius of the body. The shape factors are defined as: 


H = s*/6 (incompressible .flow) 

H = 6*/6 (compressible flow) (70) 

H x = (5 - 6*)/8 

★ 

where 6 is the boundary layer thickness, 6 the displacement thickness, 
and the emperical relationships 


H = 3.15 + - 0.01(H - l) 2 

R - 1 

2 

H + 1 5 (ti + l)|l + — ^ ^ 


(71) 


exist between H, R, and Hj. The term, r, is the temperature recovery 
factor. The entrainment coefficient is defined by 


C e = “FO— 

e e 

and the term, F, by 

(0.Q2C + C 2 + 0.8 C f /3) 

. _ e e * 0 

y = — — (0.01 vc e i 

where C- 


(72) 


(73) 


is the equilibrium skin-friction coefficient at zero pressure 
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gradient. The term C f is the skin-friction coefficient and C T is the 
shearing-stress coefficient. 

These boundary layer equations are ..integrated by a variable order, 
variable interval Adams method. The method is part of the Langley Cyber 
200 mathematical library 30 and is recomended for sets of stiff first-order 
ordinary differential equations. 

Viscous-Inviscid Interacting Theory 

Matching conditions . The global viscous-inviscid interacting 
technique depends upon a coupling of the inviscid Euler equations and the 
boundary layer equations through conventional transpiration boundary 
conditions. As pointed out by Thomas 31 , for the inviscid Euler solution 
to simulate a solution with viscous effects, it must match the vicous 
solution in that part of the flowfield where the inviscid and viscous 
equations both describe the flow accurately. The matching conditions for 
the Euler equations are described for two dimensions in his dissertation. 

The viscous-inviscid interaction technique presently used to assess 
the viscous effects on the nacelle uses a three-dimensional adaptation of 
the two-dimensional matching procedure. In both procedures, transpiration 
boundary conditions determined from a solution of the boundary layer 
equations are imposed at the body surface and in the wake of the body to 
enforce the boundary layer effects. An outline of the method used to 
determine the equivalent inviscid transpiration boundary .conditions 
necessary to match the inviscid and viscous solutions follows. 

For two-dimensional steady flow, the Navier Stokes equations are: 


and the Euler equations are 


2i+ »a. o 

3x 3y 


( 75 ) 


Designate the point where both sets of equations describe the flow 
accurately by h. Then intergrating both sets of equations with respect 
to y over the range 0 < y < h and matching the solutions at y > h 
yi el ds 



(f - F) dy 


(76) 


The subscript, w, indicates wall values and implies that the inner 
boundary of the inviscid. solution is the nacelle wall. An advantage of 
choosing the wall as the inner boundary of the inviscid solution is that 
only one computational grid needs to be generated. Thomas follows the 
example of Johnston and Sochol and lets F be a composite function 


F - f + f - f w (77) 

where f is from the boundary layer equations: 

|I + |1=0 (78) 

9x 3y 


Substituting the composite function into equation (74) and 
performing the integration gives 
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a .3 ♦ r (f - f ) dy («) 

y w y w 3x / Q w 

The resulting vector, g w , is the Euler vector, g, at the nacelle 
surface when the Euler and Navier-Stokes layer solutions are properly 
matched at y = h. Using the continuity equation as an example, the 
equivalent inviscid boundary condition at the wall necessary to match the 
inviscid and. viscous solutions is: 

(p «>» h (30) 

Note that Johnston and Sochol use the inviscid wall values in determining 
the mass fiux term on the right hand side of equation (80). This is the 
same equation as the one presented by Lock. 33 It states that the mass 
flow normal to the wall in the equivalent inviscid flow is equal to the 
streamwise rate of change of the mass flow deficit produced by the 
boundary layer. 

The previous development of the transpiration boundary conditions is 
in two dimensions. However, the present computational technique solves 
the Euler equations in three dimensions. The steady three-dimensional 
Euler equations in vector form are 


3f + 3g + 3h 
3x 3y 3z 


0 


( 81 ) 


and contain, in addition to the vector g, the vector, h, which must 
also be determined at the wall. In determining these two vectors, 
advantage was taken of the axisymmetric nature of the nacelle. By 
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performing a mass balance for -the flow between two axial stations, the 
wall, and the. edge .of the boundary layer, it can be shown that the axisym- 
metric equivalent of equation (80) is 

(PV) W {(pu) w (2r£ * + 6 * 2)} (32) 

where V is the velocity vector normal to the surface. For the three- 
dimensional adaptation, V is divided into components v_ and w_ in the 
vertical and horizontal directions respectively 


v 

-W 


U-Vw 


= d-iiA 


(83) 


Due to the relatively small boattail angles of the nacelle, the axial 
component is neglected. The vectors, n^ and n 2 , represent the y 
and z components of the unit normal to the surface. Using these 
quantities and the inviscid values of the tangent velocity, u w , and 
pressure, p w , at the surface, the equivalent inviscid vectors g w and 
h can be evaluated at the .surface of the nacelle where 

W 


pv 


pw 

puv 


puw 

pv 2 + P 

» Uw 8 ' 

pvw 

pvw 


pw 2 + p 

v(pE + p) 

w 

w(pE + p) 


(84) 
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Surface pressure equation . The. surface pressure equation for.a solid 
wall boundary with no flow across it, .equation (67), was presented 
previously in Chapter IV. For convenience, it is repeated here 



Allowing for flow across the solid boundary introduces an extra term and 
changes the surface pressure equation into 



Computations with both farms of the surface pressure equation 
resulted in negligible differences in the solution. However, all viscous- 
inviscid interacting calculations presented in this dissertation were made 
with the modified form, equation (85). 

Application of the Viscous-Inviscid 
Interaction Technique 

The transpiration boundary conditions presented in equation (84) are 
applied on both the external and internal nacelle surfaces and also in the 
wake of the nacelle. In the wake application, the transpiration boundary 
conditions are imposed along the constant C coordinate surface 
starting at the nacelle trailing edge and extending approximately one half 
of the nacelle cord downstream. At this axial location, the wake boundary 
condition had decayed to a very low value. The wake is composed of the 
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boundary layers originating from both the external -and the Internal 
nacelle surfaces. Separate transpiration boundary conditions, calculated 
from each of these components o.f the wake, were Imposed along the 
5 = constant surface depending on whether the particular grid cell was 
on the outside or inside of the surface. However, when numerically 
integrating the Euler equations,, fluxes were allowed to flow freely 
through the surface, and a pressure balance was maintained across it. 

In the overall global iteration between the Euler and boundary layer 
solutions, the boundary layer equations are solved .every 100 time steps of 
the Euler integration process using the current values of the viscous- 
inviscid solution. The transpiration boundary conditions are then updated 
using the new boundary layer solution and held constant until the next 
global iteration. The overall solution technique is started from a 
converged inviscid solution and is typically run for 1000 time steps which 
gives 10 iterations of the boundary layer. After the 10 overall global 
iterations, the solution has essentially ceased to change. 

The transpiration boundary conditions described in the preceding . 
sections physically displace the outer inviscid flow away from the surface 
to allow for the deceleration of the stream in the boundary layer. Hence 
they account for the displacement effects of the boundary layer. They do 
not account for wake curvature effects which are theoretically as 
important, but, in practice ace usually smaller in magnitude. 34 Neither 
do they account for strong interaction effects such as the interaction 
between the boundary layer and a strong shock wave which results in a 
breakdown of the usual boundary layer approxim, ' tons. 34 - 




CHAPTER VI 


INVISCID RESULTS 


Calculations were made with the two computational techniques de- 
scribed in Chapters III and IV for a flowrthrough nacelle at free-stream 
Mach numbers of 0.80 and 0.94, and at angles of attack of 0° and 4°. All 
computations were made on a Control Data Corporation Cyber 203 vector 
processor in the scalar mode. In this chapter, the results obtained with 
the alternating-direction-implicit computational procedure will be 
presented first, and then the results obtained with the explicit Runge- 
Kutta procedure. Several interesting difficulties which had a significant 
impact on the solutions were encountered and investigated in detail during 
the numerical studies. This aspect, of the research will be discussed. 

The two techniques will be compared on the basis of quality of the 
solutions, and also on practical considerations in implementing and 
processing the resulting computer codes on the Cyber 203, 


Implicit Computational Results 


Basic solution . The results obtained with the alternating-di recti on- 
implicit computational procedure are presented in figures 4 through 8. 
Figure 4 presents the basic solution for the nacelle at a free stream Mach 
number of 0.80 and an angle of attack of 0°. Wind-tunnel data of Re_and 
Peddrew^ are included In part (a) of the figure, which presents the 
Surface -pressures. A comparison of the calculated results and experiment 
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(a) Surface pressure coefficients. 


Figure 4. Basic solution calculated with the Implicit computational 
procedure. (M = 0.80, a = 0.0°. ) 
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tal data shows that the comuptational technique predicts the general 
features of the flow and the magnitude of the pressures. — 

Specifically, the calculations predict the correct axial location of 
the leading edge suction pressure peak; however, they underpredict its 
magnitude by about 20 precent. The compression and reexpansion over the 
middle portion of the nacelle is also predicted by the calculations. The 
large local gradient and reexpansion is the result of discontinuities in 
the curvature of the external nacelle surface. Although the slope of the 
nacelle surface is continuous, the juncture of the cylindrical section 
with the inlet and afterbody sections has a discontinuous curvature. 

Similar to the leading edge pressure peak, the predicted compression and 
reexpansion are in the correct axial location but the compression is some- 
what smeared by the calculations. The tendency of the computational 
technique to smear the gradients is at least partly due to the sparseness 
of the grid over the middle portion of the nacelle coupled with the 
inclusion of numerical damping in the solution technique. For example, 
over the region where the compression occurs, there are only 5 or 6 axial 
grid stations which appear to be an insufficient number to resolve the 
gradient. The sparseness of the grid may also lead to excessive numerical 

dissipation in this region. 

Part (b) of figure 4 presents computed pressure coefficient contours 
in the vertical plane for the nacelle at 0° angle of attack. The contours 
illustrate that even though the computational technique is three dimen- 
sional, the solution at 0° angle of attack exhibits the proper 
axi symmetric behavior. 
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Effect of grid refinement . A solution was obtained with twice as 
many grid points in theaxial direction as for the calculations presented 
in figure 4. The fine grid resulted in the maximum number of grid points 
possible for the implicit technique before the incore storage capacity of 
the Cyber 203 was exceeded. The solution with this, grid, which is 
presented in figure. 5, shows a considerable improvement in the agreement 
between the computations and experiment on the external surface. It has a 
more negative leading edge suction pressure peak, and the compression in 
the mid-nacelle region is stronger and has a steeper gradient. In 
addition, the computed reexpansion region on the nacelle is in quite good 
agreement with the wind-tunnel data.. 

Nacelle at angle of attack. A solution for the flow around the 
nacelle was also obtained at a free -stream Mach number of 0.80 and an 
angle of attack of 4.0°'. Figure 6 presents a comparison of the calculated 
pressures with the experimental data of Re and Peddrew for the side 
meridian of the nacelle, <f> = 90°. The computational technique predicts 
the qualitative and quantitative character of the flow well, and the 
general comments about the comparison of the calculations with 
experimental data at 0° angle of attack apply. For this calculation, 
however, it was necessary to average the internal and external flow 
quantities at the leading edge in order to obtain a solution. 

The computed pressure distributions for the top, side, and bottom 
rows of the nacelle are presented in figure 7. The external distribu- 
tions, presented in figure 7(a), demonstrate the threer^ijmensional 
character of the flow with the nacelle at angle of attack. The majority 
of the three-dimensional effect is confined to the forward portion of the 
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(a) Outside surface. 

Figure 7. Computed surface pressure coefficient distributions for 
several rows at angle of attack. (M w * 0.80, a « 4.0°, 
implicit code.) 
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nacelle where the calculations show that a large suction peak is being 
generated at the leadi.ng edge of the top meridian. The internal pressure 

coefficients, presented in figure_7(b), are influenced by the three*?.. 

dimensionality of the flow only in a relatively small region at the 
nacelle leading edge.. 

Contours of the computed pressure coefficients in the vertical plane 
with the nacelle at an angle of attack of 4.0° are presented in 
figure 8. The contours also illustrate the three-dimensionality of the 
computed flowfleld with the nacelle at angle of attack. A comparison of 
these contours with those at 0.0° angle of attack, presented in 
figure 4(b), illustrates the more pronounced pressure gradient on the 
inside of the nacelle at a = 0.0°, and the differences, between the 
pressure gradients on the inner and outer surfaces near the leading 
edge. It also illustrates the greater expansion of the external flow on 
the top of the nacelle at the higher angle of attack. 

Explicit Computational Results 

Solutions for the flow past the flow-through nacelle were also 
computed using the explicit Runge-Kutta computational procedure with 
implicit smoothing of the residuals described in Chapter IV. The results 
obtained with the explicit procedure are presented in figures 9 through 
12 . 


Basic solution . The basic solution at a. Mach number of 0.80 and an 
angle of attack of 0° is presented in figure. 9. Like the implicit 
computational technique, the explicit technique predicts all of the 
general features of the flow, and the pressures on the external surface 
agree quite well with the wind-tunnel data (see part (a) of the figure). 



Figure 8. Pressure coefficient contours for the computations at angle 

of attack. (M. = 0.80, a = 4.0°, vertical plane, implicit code.) 

op 
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(a) Pressure distributions. 


Figure 9. Basic solution computed with the explicit computational 
procedure. (M = 0.80, a = 0.0°. ) 
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By comparing these pressure distributions with the Ones calculated by 
using the implicit procedure, figures 4(a) and 5, one can see that the 
explicit procedure gives a more accurate solution on the external surface 
for the same number of grid points. This is particularly evident in the 
region of the rapid compression and reexpansion which occurs on this 
surface. Here the explicit procedure captures the gradients much better 
than the implicit procedure and comes close to matching the resolution of 
the pressure distribution computed by using the implicit procedure with 
twice as many grid points. The computed pressures on the internal surface 
are much more positive than the data. 

Pressure .coefficient contours for the explicit solution are 
presented in part (b) of figure 9. The region of low pressures which is 
affected by the nacelle is slightly larger for the explicit calculations 
than for the implicit calculations (see figure 4). The differences in 
both the external surface pressures and the pressure coefficient contours 
between the two solutions may be the result of the slightly different 
implementation of the surface boundary condition. 

Grid refinement study . An attempt was made to look at the effect of 
grid resolution for the explicit computational technique. It was possible 
to add 10 grid stations, axially along the surface of. the nacelle, for a 
total of 40 grid stations along the surface, before the incore storage 
capacity of the Cyber 203 was exceeded. Figure 10 presents the solution 
for this finer mesh at a free stream Mach number of 0.80 and an angle of 
attack of 0.0°. By comparing the computed pressures, which are presented 
in figure 10(a), with the pressures for the coarser mesh, one can detect a 
slight improvement in the agreement with the wind-tunnel data as the mesh 
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(a) Pressure distributions. 

Figure 10. Solution for the finer- grid and the explicit computational 
procedure. (M = 0..80, a = 0.0°. ) 
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Is refined, the pressure coefficient contours for the refined grid case 
are presented in part (b) of figure 10, ... The contours axe. essentially the 
same a$ those for the coarser grid. 

The modest increase in. the number of grid points resulted in a 
slight improvement in the correlation between the explicit solution. and 
the wind-tunnel data. Even though the overall level of the correlation is 
very, good, a greater increase in the number of grid points may further 
improve the agreement in the region of the peak expansion and the 
following compression. 

Internal pressure level . The addition of more grid points along the 
surface did little to improve the correlation between the computations and 
the wind-tunnel data on the internal surface. The results of the implicit 
technique agree much better with experiment inside the nacelle, and hence 
seem to be more accurate. However, the implicit technique has a Kutta- 
like condition imposed at the trailing edge. A calculation with a version 
of the implicit code which did not contain the "Kutta" condition shows the 
same basic, trends as the explicit code as can be seen from figure 11. In 
a similar observation, Miranda^ indicates that inviscid potential 
solutions for similar configurations frequently give good solutions for 
the external flow but yield the incorrect internal mass flow ratio. These 
results lead one to speculate that the Euler "no.n-Kutta" solutions possess 
the correct internal trends for inviscicLflow. Chapter VII examines in 
detail the. possibility that the discrepancy between the data and the 
calculations is due to viscous dissipation present In the experiment. . 

Supercritical solution . In addition to the solution at a Mach number 
of 0.80 and an angle of attack of 0°, a solution with the explicit 
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Figure 11. Comparison of experimental and calculated p-essures for the 
inviscid Euler equations. (M = 0.80. a = 0°.l 







procedure was obtained- at a Mach number of 0.94 and 0°..angle of attack. 

The experimental data -indicated that 0.94 was the highest Mach number for 
which the flow on the nacelle was not separated. The calculation was made 
with, the finer computational mesh described in the earlier, section on grid 
resolution . 

The computational results agree very well with the measured pres- 
sures on the external surface of the nacelle as figure 12(a) 
illustrates.. On the internal surface, just as they did at the Mach number 
of 0.80, the calculations predict higher pressures than the wind-tunnel 
data. Even though the predicted peak expansion near the leading edge is 
low, the pressures agree with the data everywhere except at the very 
peak. Probably grid resolution or excessive dissipation due to the large 
changes in the pressures in this region of the flow is responsible for the 
discrepency. The calculations both qualitatively and quantatively predict 
the rearward movement of the leading edge negative pressure peak with 
increasing Mach number. Notice that the general shape of the leading edge 
expansion and subsequent recompression has changed from the lower Mach 
number. At the higher Mach number, the expansion continues until it is 
abruptly terminated by a strong compression, or possibly a shock. The 
calculations correctly reflect this shape change. Also, at the higher. 
Mach number, the more pronounced hump in the pressures on the rear of the. 
nacelle is predicted and the calculations match the data very well in this 
regi on . 

Parts, (b) and.(c) of figure 12 present pressure contours in the 
vertical plane for the calculation. Part (b) of the figure shows the 
overall region in the vicinity of the nacelle, and -part (c) presents a 
detailed view of the leading edge region. The pressure contour 
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(a) Pressure distributions. 


Figure 12. Supercritical solution calculated with the explicit 
computational procedure. (fl = 0.94, a = 0.0°. ) 
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corresponding to a local Mach number of 1.0 is represented by a dashed 
line. The. contours illustrate the extent of the supersonic bubble On the 
forward part of the nacelle, and also show that there is a small 
supersonic bubble on the afterbody. In addition, they, illustrate that the 
zone influenced by the nacelle, both the compression zones originating at 
the leading and trailing edges and the expansion zones, extend further 
outward in the radial direction. One might expect this general change in 
the nature of the flow at the supercritical Mach number. 

The calculations have correctly predicted the changing nature of the 
flow from the lower to the higher transonic Mach number, and also 
predicted the quantitive results at the two Mach. numbers very well. 
Combined with the good predictions at angle of attack, these results 
demonstrate the potential of the Euler equations in solving flows of this 
complexity. 

Numerical Problem Areas 

During the development of the implicit computational procedure, and 
the subsequent numerical studies using both the implicit and explicit 
procedures, two unexpected numerical difficulties were discovered and 
investigated. The first difficulty concerns the stability properties of 
the implicit algorithm. The second pertains to a numerically produced 
surface total pressure loss, and -is inherent in both the implicit and 
explicit algorithms. While other researchers may have encountered similar 
difficulties, until recently a large segment of the computational 
community was unaware of the .stability problem. 37 The total pressure 
problem still remains largely un reported. 37 The Investigation of these 
areas is discussed in this section. 
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Effect of boundary conditions on implicit convergence . In the early 
ISBO's, when the development of the Implicit computational procedure was 
initiated, .the Beam and Warming alternating-direction-implicit numerical 
algorithm was relatively new and untested, particularly in three 
dimensions. Difficulties with its stability properties in three 
dimensions were then known or. suspected by only a small community of 
researchers. 38 " 39 As noted in Chapter III, while .the implicit phase of 
the present work was eing conducted, Abarbanel , Dwoyer, and Gottleib 23 
proved that the undamped Beam and Warming scheme in three dimensions is 
weakly, but unconditionally, unstable. They also showed that the weak 
instability is only present in the very long wavelengths. 

During the development of the implicit procedure, an instability 
manifested itself near the outflow boundary. The problem was 
investigated, and a radiation outflow boundary treatment was applied which 
enabled reasonably accurate engineering solutions to be obtained. These 
results which are reported by Compton and Whitesides in reference 14 tend 
to confirm the work of Abarbanel, Dwoyer, and Gottleib. 23 The results of 
reference Is have .helped clarify the weak instability in the undamped Beam 
and Warming algorithm when applied, to the three dimensional Euler 
equations. Thie aspect of the development of the implicit computational 
procedure is described below. 

Two parameters, AC p>max and the maximum residual , were used to 
test the convergence of the numerical solutions obtained with the Implicit 
numerical technique. At any given time level, ACp Is defined to be the. 
absolute change in the pressure coefficient between the present time level 
and the previous time level. The total residual at any time step is 
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defined to be the root mean square of -the rate of change of all five of 
the normalized dependent variables. At. every time step, ACp and the 
total residual are computed at each grid point and the entire 
computational region is searched for the maximum value of each parameter. 

Figure 13 presents the iteration history for the solution obtained 
with the outflow boundary condition p = p # and extrapolation at the far- 
field boundary. The history of AC indicates that the solution is 

p f ilia a 

converging. However, the history of the maximum residual indicates that 
the solution is in fact diverging. Thus it is misleading to base 
convergence strictly o_n the change of a flow parameter. The increase in 
maximum residual coincides with its location gradually changing from the 
vicinity of the nacelle to the outflow boundary, implying that the outflow 
boundary condition p = p w equation (25), is partially reflective. 

The nonceflecting outflow boundary condition of Rudy‘S was tried 
with mixed results in stabilizing the solution process. In addition, the 
Riemann-invari ant method of treating the outflow boundary,^ which was 
used very successfully in combination with the explicit computational 
procedure, was briefly investigated. However, a satisfactory solution was 
not obtained with this combination. 

As a consequence .of these results, the radiation boundary condition, 
equation (32), was derived and investigated. The iteration history of the 
residuals when the radiation boundary condition was used at both the 
outflow and far-field boundaries is presented in figure 13(b). The 
overall maximum residual Indicates that the pressure disturbances passed 
through the outflow boundary more easily, delaying the emergence of the 
instability until approximately 1000 iterations. After the instability 



ITERATIONS ITERATIONS 

(a) Outflow boundary condition p = p^. 

Figure 13. Iteration history for the residual of the implicit 
computational procedure. (M = 0.80, a s 0.0°.) 
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(b) Radiation outflow boundary condition. 


Figure 13. Concluded. 
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emerges, the rate at which the overall maximum, residual grows Is reduced 
to near zero from the rapid growth associated with using equation 25, 

Validity of the implicit solutions . The Iteration history of the 
overall maximum residual presented in figure 13(b) tends to confirm the 
results of Abarbanel et. al . However, the boundary conditions were 
imposed explicitly which may degrade the convergence properties of the 
implicit scheme and increase the run time. The maximum residual near the 
nacelle is also plotted in figure 13(b). This local residual continues to 
decrease with a total drop of about 4 orders of magnitude in 1200 
iterations, and 6 orders of magnitude in 2400 iterations, further 
indicating that the instability is associated with the boundary treatment. 

Residual contours after 1200 iterations or time steps are 
presented in figure 14. These contours confirm that the maximum residuals 
do indeed occur near the outflow boundary and that the solution near the 
nacelle appears to be converging. Since the solution near the nacelle 
continues to converge and the local residual has decreased 4 orders of 
magnitude at 1200 time steps, the solutions near the nacelle should be 
reasonably accurate and useful for engineering calculations. 

Explicit convergence properties . There does not appear to be a 
problem with the stability or convergence of the explicit computational 
procedure. The presentation of these properties of the explicit technique 
is placed in this sectior merely for comparison with the implicit 
algorithm. The residual upon which convergence is based for the explicit 
computational procedure is the difference between the computed enthalpy at 
any particular grid point and the free stream enthalpy. Figure 15 
presents the iteration .history Of the maximum residual in the 
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Figure 15. Iteration history for the residual of the explicit 
computational procedure. (M w = 0.80, a = 0.0°. ) 
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computational region. for the bas.ij: solution. The history shows that the 
solution is stable and converging. . The maximum residual Initially drops 
rapidly,, and altnough its rate of decrease levels off some, its general 
trend_still continues down after 4000 time steps. 

The solutions were obtained at a Courant number of 4 which was the 
largest Courant number for which numerical stability could be maintained. 
The optimum rate of convergence may be better than indicated by 
figure 15. However, no studies were made to determine the size of the 
damping parameters for the maximum rate of convergence. Numerical studies, 
by Vatsa 40 indicate .that the optimum values of the damping parameters for 
a maximum convergence rate depends upon the grid. 

Total pressure loss at the surface. Cumulations with both computer 
codes show a surface total pressure loss even at subsonic Mach numbers. 

This feature of the solutions, which is inconsistent with the physics of 
the inviscid flow, was first noticed in calculations made with the 
implicit computer code, and was reported in reference 14. However., an 
indepth investigation of the problem was not. undertaken at the time. When 
solutions obtained with the explicit code also exhibited this 
characteristic, it was considered highly desirable to investigate it more 
fully. The resulting investigation of the problem is discussed in this 
section. 

Since there are no supersonic regions, and hence no shocks in the 
present .solutions for * 0.80, and since the flow is Inviscid, the 
total pressure in the entire flow field should be constant. However, the 
numerical calculations at this subsonic Mach number result in total 
pressure losses which are most noticeable on the Internal surface of the 
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nacelle. Although the losses are. confined to a small region in the 
immediate vicinity of the..nacelie, they, are sometimes quite large. The 
total pressure losses do not seern to affect the values of the static 
pressures, but as stated before, they are inconsistent with the correct 
physical nature of the problem. _ 

Figure 16 presents the total pressure distributions on the nacelle 
surface, and illustrates the magnitude of the total pressure losses for 
the free stream Mach number of 0.80. Tne figure is typical for both the 
sparse and fine grids since refining the grid-had little effect on the 
losses. Parts (a) and (b) of figure 16 present the surface total 
pressures calculated by using the implicit numerical procedure. The 
results presented in figure 16(a) were obtained with an. early version of 
the code in which the flow variables for the top and bottom surface were 
averaged at theJeariing edge. While. the total pressure on the external 
surface of the nacelle is essentially correct, the figure shows a loss in 
total pressure of approximately 5,0 percent on the internal surface of the 
nacelle. Initially, the loss was attributed to incorrect treatment of the 
the leading edge boundary (the implicit computational technique required a 
direct treatment of the leading and trailing edge boundaries). As a 
consequence, a wide variety of leading edge treatments was Investigated in 
order to correct the problem. 

The outcome of the leading edge study was to adopt the present 
treatment described in Chapter III in which the two surfaces were treated 
separately. Figure 16(b) shows the total pressure distributions resulting 
from the calculations using the new leading edge boundary condition. The 
loss Is now split between the upper and lower surfaces, with a loss of 
about 2.5 percent On each surface. Considering the following results 
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(a) Implicit computitional procedure, dependent variables 
averaged at the leading edge. 

Figure 16. Surface total pressure distributions. 

(M = 0.80, a * 0.0°.) 
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(b) Implicit computational procedure, dependent variables not 
averaged at the leading edge. 


Figure 16. Continued 





obtained with the explicit computational procedure, it Is doubtful that 
the new boundary treatment, results in a more -physically correct solution. 

Since the explicit computational procedure incorporates a finite 
volume differencing technique, and therefore does not require any direct 
treatment of the leading or trailing edge boundaries, it was not expected 
to give total pressure losses at a subsonic Mach, number. However, as part 
(c) of figure 16 shows, the problem is magnified by the explicit 
technique. Like the results for the implicit technique with the earlier 
leading edge boundary treatment, the total pressure on the external 
surface is essentially correct. The losses on the internal surface, 
however, have increased to approximately 10 percent. These results 
indicate that the source of the loss is not associated with the treatment, 
of the leading and trailing edge boundary conditions. 

Origin of the total pressue loss . In order to help determine the 
source of the losses, additional numerical experiments were conducted with 
the explicit computational procedure using an unswept, untapered wing as 
the test configuration. Using the unswept, untapered wing resulted in a 
physical problem with two-dimensional flow, but retained the three- 
dimensionality of the computational procedure. .The results of the 
experiments, which are presented in the appendix, indicate that the 
surface boundary condition as well as the rest of the explicit computer 
code was programmed correctly. The numerical studies on the wing also 
indicate that the origin of the total pressure loss Is. an expansion of the 
flow around the sharp leading edge of the nacelle. The sharp leading edge 
is a feature of the nacelle geometry made necessary by the H-grid. The 
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studies additionally indicate that by eliminating the nearly discontinuous 
expansion, the total pressure losses can be .eJ iminated. 

These ideas are demonstrated for the nacelle in figures 17 and 18. 
Figure 17 presents the solution for the computations made with the 
explicit numerical computational technique at a free stream Mach number 
of 0,80 and an angle of attack of 0.0°. Parts (a) and (b) of the figure, 
show respectively the Mach number contours and the velocity vectors in the 
vicinity of the leading edge of the nacelle. They illustrate that 
contrary to what one would expect at 0.0° angle of attack, the flow 
stagnates on the external surface of the nacelle and expands around the 
sharp leading edge to the internal surface. 

The stagnation point on the external surface can be seen most 
clearly from the Mach number contours. The clustering of the Mach number 
contours at the leading edge indicates the severity of the singularity 
created in the flow as it expands in traversing the discontinuity in the 
nacelle surface. The velocity .vectors show an inward component of the 
flow immediately above and in front of the leading edge, and futher 
indicate that the direction of the flow at the leading edge is from the 
external surface, around the leading edge, and to the internal surface. 
Since the first-order damping terms become very large .in regions with such 
strong gradients in the flow,, a calculation was made. with these 
dissipation terms turned off. Very little difference could be detected in 
the solution. 

It is .thought that this problem would be much less severe for a 
rounded leading edge coupled with the use of a locally embedded body 
normal grid such as those normally used for blunt- bodies . However, as the 
appendix shows, care must be taken even with a rounded leading edge and a 



(a) Mach number contours. 


Figure 17 


. Calculated flow-field in the vicinity of the 
leading edge of the flow-through nacelle. 

(M = 0.80, c» * 0.0°, explicit procedure.) 





C.-type grid. Even with, this type grid,. the. strong gradients generated as 
the flow expands around the rounded leading edge can for all practical 
purposes appear as dlscontlnultes unless there are sufficient grid points 
to resolve the rapid changes in the flow variables. ........ 

Eliminating the total pressure loss . To demonstrate that the total 
pressure losses on the internal surface of the nacelle can be eliminated 
by eliminating the expansion around the sharp leading edge, calculations 
were made for a modified flow-through nacelle. The modification to the 
nacelle consisted of adjusting the mean camber line of the nacelle's air- 
foil so that the flow stagnated precisely at the leading edge instead of 
on the external surface.. Thus the expansion of the flow around the sharp 
leading edge was eliminated. 

The camber was decreased to 75 percent of the original camber at the 
leading edge. Along the airfoil, the reduction in camber was 
progrpssi vely decreased according to the square of the cord so that at the 
trailing edge, the camber remained 100 percent of the original value. The 
change in the mean camber line necessary to.adjust the stagnation point 
was most noticeable in the first 25 percent of the cord. Figure 18 
presents the solution for the modified nacelle with part (a) of the figure 
showing the pressure coefficient (List ri but ions. The pressures look 
reasonable; however, comparisons with experiment. could not be made since 
•no wind-tunnel data was available for the modified configuration. 

Mach number contours in the vicinity of the leading edge are 
presented in part (b) of the figure, they illustrate the fact that the 
flow dues stagnate precisely at the leading edge, and hence does not 
expand around around. the corner. The result, shown in figure 18(c), is 



(a) Surface pressure coefficient distributions. 


Figure 18. Solution for... the modified flow-through nacelle 
(M w 0.80, a = 0.0°, explicit procedure.) 
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(c) Surface total pressure distributions 


Figure 18. Concluded 


that the total pressure on both the external and the internal surfaces of 
the nacelle maintains Its free stream value. An Indication of the amount 
of change in the nacelle's mean camber line necessary to adjust the 
stagnation point can be seen by comparing the modified nacelle contour 
depicted in figure 18(b) with the original nacelle contour depicted in 
figure 17(a). 

The source o' the total pressure loss on the surface of the nacelle 
has been identified as an expansion around the discontinuity in the 
nacelle's surface at the .leading edge. It has also been demonstrated that 
by eliminating the expansion around the sharp leading edge, the free 
stream total pressure can be maintained, . 

Comparison of Techniques 

The previous sections in this Chapter have presented a discussion of 
the solutions calculated by the two computational techniques, an implicit 
technique employing the three-dimensional Beam and Warming numerical 
algorithm, and an explicit technique employing the four-stage Runge-Kutta 
algorithm with implicit residual smoothing. The quality of the .solutions 
obtained with these techniques was compared during the discussion.. The 
comparison on the basis of solution quality will be summarized in this 
section, and, in addition, practical considerations of implementing the 
two techniques on the computer will be discussed. 

Processi ng time and computer storage requi rements . All numerical 
calculations were processed on a Control Data Corporation Cyber 203 
computer in the scalar mode. For the implicit code, a typical calculation 
was processed for approximately 1100 time steps on a computational grid 
containing 18,502 grid points, the calculations were made at a Courant 
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number of Orde£ 11 based on the minimum grid spacing in the axial 
direction. .The processing time on the Cyber 203 was approximately 6.7- 

hours or 1.22 milliseconds per grid point per time step.. The code was .... 

written for maximum visibility of the form of the equations, and hence 
contained many divisions by quantities whose value, were 1.00. By 
eliminating these divisions, and optimising the code for speed of 
computation, a considerable reduction in processing time should be 
real i zed. 

The explicit code required an average of 3000 time steps to reach a 
steady state solution. For the same computational grid, this yielded a 
computer processing time of approximately 2.1 hours or 0.16 milliseconds 
per grid point per. time step or approximately one third the computational 
time of the implicit code. Although neither -the implicit code nor. the 
explicit code was vectorized, vectorization of the codes should result in 
a decrease in the processing time. 

In order to calculate flows around more complicated configurations, 
additional grid points will be required. The maximum number of grid 
points possible for the implicit computer code as it was written was 
approximately 37,000 before the incore storage capacity of the Cyber 203 
computer was exceeded. By. deleting all possible arrays in the computer 
program, and by recomputing the necessary parameters from the deleted 
arrays every time they were needed, approximately 63,000 grid points could 

be obtained before the computer's incore storage capacity was exceeded. 

The maximum number of grid points possible before. the explicit code, which 
was not_optim1zed for minimum storage requirements, exceeded the incore 
storage capacity of the CyDer 203 was approximately 20,000. 


In general, the Implicit code requires less computer storage while 


the explicit code runs faster. Both the run time and the computer storage 

requirements of each code are marginal for general design studies where 

many configurations need to be evaluated quickly unless optimization of 
programming for storage and run time can be realized. However, the cost 
of obtaining a given, solution has steadily decreased by an order of 
magnitude every eight years for the past fifteen years. 41 Computers which 
have much more computing power and memory than the Cyber 203 will bebome 
available in the very near future. For example, NASA Langley Research 
Center .recei ved a 16 million word version of the Control Data Corporation 
Cyber 205 in the fall of 1984. Another example is the MAS facility 
discussed by Ballhaus. 41 With such advances in computer technology, 
coupled with advances in numerical techniques, computer codes based on 
solutions of the three-dimensional Euler equations should soon become 
practical for problems of this complexity. 

Convergence . While reasonably accurate engineering solutions were 
obtained with the implicit technique, its overall convergence properties 
are unacceptable. The implicit application of the boundary conditions to 
this problem may enhance its stability and should be investigated, as 
should the effect of enthalpy damping on its stability, characteristics. 
However, at present, the scheme has a weak instability when applied to 
this three-dimensonal problem with the result that the implicit solution 
never actually converges. The explicit scheme Is stable and converges 
even though Its rate of convergence is less than desirable. .. 

Accuracy . The relative accuracy of the solutions obtained with the 
two numerical techniques has already been discussed in the previous 
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•sections. Basically, the explicit technique- yields more accurate 
solutions for the same number of grid points in that it captures the 
pressure gradients and peak pressure expansions better.than the implicit - 
technique. The better accuracy of the explicit code may be a result of 
the slightly different implementation of the surface boundary conditions 
in the two numerical procedures. Both procedures, however, give solutions 
that. are in good agreement with experimental data on. the external surface 
of the nacelle where the effects o f viscosity are relatively small. 



CHAPTER VII 


VISCOUS-INVISCID INTERACTING RESULTS 

In. general, the inviscid solutions calculated by using the Euler 
computational techniques agree very well with the wind-tunnel data on the 
external surface of the nacelle. However, there are large differences 
between the inviscid solutions and the measured pressures inside the 
nacelle's duct.. In Chapter VI, it was hypothesized that since the Euler 
equations model compressibility and rotational ity, the discrepancy was due 
to viscous dissipation present in the experiement. In the physically 
realistic case, strong interactions often occur between the viscous 
boundary layer and the main. stream even when the boundary layer does not 
separate. Typical examples of this type of flow occur near the trailing 
edge of an aft-loaded airfoil or at a shock -boundary-layer interaction. 

In these instances, modeling the frictional forces becomes essential if 
accuracy is to be maintained. 

Therefore, to complete the study of the solution of the flow field 
for the flow-through nacelle, a preliminary assessment was made of the the 
influence of ' viscosity. The viscous-inviscid interacting computational, 
model described in Chapter V was used for this phase of the research which 
was conducted at a free stream Mach number of 0.80 and an angle*of-attack 
of 0.0°. These investigations, which have yielded new Insight into the 
mechanics of the interactions between the Internal and external flows, are 
described in the present chapter. 
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Results 


Overall viscous effects . The overall result of Including the 
boundary layer effects in the otherwise inviscid Euler solution is illus- 
trated in figure 19. The figure presents the viscous-inviscid interacting 
solution as well as the inviscid solution, and illustrates that the 
viscous effects significantly change the internal pressures but leave the 
external pressures largely unchanged.. Including the boundary layer and 
wake does. result in a slight decrease in the external pressures very . near 
the trailing edge.... However, on the internal surface, it causes a 
significant decrease in the exit pressure, and produces a sizable axial 
pressure gradient in the nacelle's duct. The net effect improves the 
correlation between the computed internal pressures and the experimental 
data. 

As pointed out previously, in addition to surface boundary layer 
thickness, this interacting procedure compensates for wake thickness but 
not wake curvature. Melnik^ indicates that wake curvature,, while not 
being as important as wake thickness, produces similar results. Hence, if 
allowance were made for wake curvature in the present calculations, the 
computations should match the internal, pressure datum even better. 

External and internal boundary layers . An attempt was made to 
determine the relative Importance of the external and internal boundary 
layers in changing the character of the flow past the nacelle. Therefore, 
in addition to the Inviscid and the complete viscous-inviscid interacting 
solutions, interacting computations were made with the the boundary layer 


and wake originating only from the external surface of the nacelle. 
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Figure 19. Effect of viscous-inviscid interaction on the 
nacelle pressures. (M * 0.80, a = 0.0°. ) 
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Similar commutations were also made with the boundary layer and wake 
originating only from the Internal surface. 

Figure 20 presents the. computed pressure distributions on the 
external surface of the nacelle for all four solutions, and shows that 
neither the external nor the internal boundary layers greatly influence 
the '..tsiJe pressures. .The viscous effects on the internal nacelle pres- 
sures are illustrated in figure 21. These effects are substantial, and 
consist of both a change in the overall pressure level and the generation 
of a pressure gradient in the axial direction. 

First consider the influence of the boundary layer and wake 
originating from. only the external surface of the nacelle. Part (a) of 
figure 21 contains the pressure distributions for the complete viscous- 
inviscid interacting solution, and for the solution which includes the 
viscous effects produced by only the external nacelle surface. The common 
factor between the curves is that each includes the viscous effects 
originating from the external surface of the nacelle. Even though there 
are differences in the pressure gradients between the two calculations, 
the exit pressures in both cases are the same. 

In part (b) of figure 21, the viscous effects due to the external 
nacelle surface are absent. Part.(b) presents the inviscid solution and 
the interacting solution with only the internal boundary Jayer and wake. 

As in part (a), the pressure gradients are different but the exit, 
pressures of both curves are equal. The message of the ..comparisons is 
that the external boundary layer and its wake in combination with the 
inviscid flow determine the exit pressure and hence the overall pressure 
level of the internal flow. The internal boundary layer has very little 
effect on the exit pressure. Comparing parts (a) and (b) also illustrates 
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(a) Total versus external viscous-inviscid interaction. 



(b) Inviscid versus internal viscous-inviscid interaction. 


Figure 21. Viscous effects on the internal nacelle pressures. 
(M* = 0.80, a = 0.0°. ) 
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(c) Inviscid versus external viscous-inviscid interaction. 
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(d) Total versus internal viscous-inviscid Interaction. 


Figure 21. Concluded 
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thiit the external viscous effects produce a lower exit pressure and hence 

a lower overall internal pressure that Is in better agreement wi t h the 

wind-tunnel data. 

Next, consider the effect of the boundary layer originating from 
only the internal surface of the nacelle. In. parts (c) and (d) of . 
figure 21, the correlating factor, is the viscous effects produced by the 
inside nacelle surface. In part (c), both pressure distributions are the 
result of calculations in which the internal viscous effects are absent, 
and in part, .(d) both are the result of calculations with them present. An 
examination of both sets of pressure distributions illustrates that the 
boundary layer on. the internal surface of the nacelle produces a pressure 
gradient in the nacelle duct. In addition, it shows that the gradient is 
essentially unaffected by the boundary layer on the outside surface. One 
dimensional axisymmetric calculations demonstrate that the gradient is the 
result of the change in the effective duct area due to the growth of the 
internal boundary layer. For example, the one-dimensional calculation 
yielded a pressure coefficient gradient of 0.13 versus the gradient of 
0.11 given by the present viscous-inviscid calculation. 

Interacting mechanism . A more complete understanding of the 
mechanism by which the boundary layer produces these results, and of the 
relative influence of each boundary layer..on the overall flow pattern can 
be gained by also considering figures 22, 23 and 24. Figure 22 shows the 
strength of the sources and sinks that are Imposed on the nacelle surface 
and in the wake in. order to satisfy the transpiration boundary conditions 
for the continuity equation. The strengths are equivalent to the stream- 
wise rate of change -of the mass flow deficit produced by the boundary 
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layer or wake. The figure Illustrates that the external boundary layer, 
and Its associated wake result in boundary conditions that are nearly an 
orde'' of magnitude larger than those associated with the internal boundary 
layer and wake. The figure further shows that the dominating portion of 
the external viscous effects comes from the trailing edge region where 
there is a rapid deceleration of the external flow and consequently a 
rapid buildup of the boundary layer. 

The influence of these viscous boundary conditions on the flow field 
in the vicinity of the trailing edge is illustrated in figure 23. Part 
(a) of the figure shows the velocity vectors in the vicinity of the 
trailing edge for the inviscid solution, and part (b) the corresponding 
vectors for the interacting solution. Comparing the two velocity vector 
plots illustrates the difference in the basic nature of the two solutions; 
the inviscid solution possesses a greater inward radial component of the 
flow which suggests a greater circulation. 

A quantitative comparison of the viscous boundary conditions on the 
velocity vectors. immediately downstream of the trailing edge of the 
nacelle is presented in figure 24. Part (a) of the figure presents the 
velocity vectors of the inviscid Euler solution, and the ' interacting solu- 
tion with only the external viscous effects. It shows that che viscous 
effects originating. from the external surface straighten out the. 
streamlines, thus reducing the circulation in the trailing edge region. 

The effects of viscosity produced by the internal surface, part (b), show 
the opposite effect, but are much smaller in magnitude as would be 
expected from the results presented in figure 22. The net effect on the 
velocity vectors immediately downstream of the trailing edge is presented 
in figure 24(c). Basically, the viscous effects reduce the inward radial 



(b) Inviscid versus internal viscous-inviscid interaction 


Figure 24. Detailed comparison of velocities immediately 
downstream of. the trailing edge, (viscous and 
inviscid solutions, M = 0.80, u = 0.0°. ) 
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velocity component of the mainstream and increase its streamwise 
component . 

Decambering concept . The external viscous effects produce a change 

in the stream velocities that is similar to a change that would be 
produced by decambering the nacelle's airfoil. The analysis presented 
above strongly indicates that the decambering effect is responsible for 
decreasing the exit pressure and the overall internal pressure level from 
their inviscid values. This theory can be further evaluated by two 
completely inviscid. tests. For the first test, recall the comparison of 
the two inviscid solutions which were calculated with the implicit 
computational procedure (see figures 3(a) and 11). For. the "non-Kutta" 
calculation presented in figure 11, a substantial inward radial velocity 
existed at the trailing edge of the nacelle. For the "Kutta" calculation 
presented in figure 3(a), this velocity was set equal to zero. As would 
be expected, imposing the "Kutta" condition . lowered the internal pressure 
level so that it agrees much better with experiment. 

The second test was conducted with the explicit Euler code. In this 
test, an infinitely thin trailing edge extension, or tab, was attached to 
the nacelle parallel to its axis of symmetry. Calculations using this 
configuration produced even more dramatic results than those of the first 
test. The -internal pressure coefficient was lowered from 0.3, a value 
which was higher than the experimental pressure coefficient, to 0,1, a 
value approximately one half the experimental pressure coefficient. 


The results presented above show that the external boundary layer 
effectively deCambers the airfoil of the nacelle and alters the overall 
flow pattern by redirecting the streamlines closer to the free stream 
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direction, A sketch of.tne inviscid and viscous streamlines is presented in 
figure 25., The compression on the external surface of the nacelle is. thus 
weakened, resulting in a less positive exit pressure and a less positive 
internal pressure level which agrees better with wind-tunnel data than the 
iv/iscid computations. 

Implications of viscous effects . It is obvious that the. magnitude 
of these results is configuration dependent. However, one concludes that, 
in. the absence of any artificial "Kutta" condition, the Euler equations 
give the correct i nviscid .trends . That is, that the i n viscid pressure 
level inside the nacelle's duct is more positive than would be experienced 
in reality. Hence, if simulation of the correct mass flow through the 
nacelle's duct is important in analyzing, a fluid flow problem, then 
viscous effects must be included, in the computational model. 

These implications also extend to wind-tunnel testing. of models with 
flow-through nacelles. In the early stages of developing a new aircraft, 
when many different configurations must be tested,, the complications and 
expense associated with. testing powered models necessitates the use of 
models with flow-through nacelles. Adjustments to the flow-through data 
are tnen determined by testing powered versions of the most promising 
configurations. This experimental procedure was followed by Capone 
et. al .? in determining the overall, aerodynamic characteristics of a new 
fighter aircraft. As the present viscous-invisctd investigation shows, 
the external boundary layer on the nacelles can affect the internal mass 
flow as well as the pattern of the streamlines both behind and in front of 
the nacelle. Thus, the accuracy of simulating the nacelle boundary layer 
could significantly impact the measured aerodynamic characteristics of the 
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configuration. 



Summary of the Viscous Effects 

In summary, for predicting the flow on the external surface of the 
nacelle, viscous effects were relatively unimportant and the Euler 
equations gave good solutions. In contrast, for predicting the flow 
inside the nacelle's duct, the viscous effects are extremely important and 
both the external and internal boundary layers and wakes must be 
simulated. The internal boundary .layer creates an axial pressure gradient 
in the nacelle's duct, but essentially does not affect the overall 
pressure level. The external bound.' r y Jayer and its associated wake 
change the overall pattern of the inviscid flow. They displace the 
streamlines away from the external surface of the nacelle thus effectively 
clambering the nacelle's airfoil. As a result, the compression at the 
trailing, edge is weakened. This gives a less positive exit pressure, and 
hence a less positive overall internal pressure level which is closer to 
the free stream value and agrees better with wind-tunnel data than the 
inviscid computations. Hence^ if simulating the correct mass flow through 
the nacelle's duct is important, then viscous effects must be included in 


the computational model. 



CHAPTER Vi 1 1 


CONCLUDING REMARKS 

A study has been made of the solution of the three-dimensional flow 
field for a flow-through nacelle. Both inviscid and viscous-inviscid 
interacting solutions were examined. Inviscid solutions were obtained 
with two different computational procedures for solving the three- 
dimensional Euler equations. The first procedure employs an 
approximately-factored alternating-direction-implicit numerical algorithm, 
and required the development of a complete computational model 
specifically geared to the nacelle problem. The second computational 
technique employs a fourth -order Runge.-Kut.ta numerical algorithm which was 
modified . to fit the nacelle problem. Viscous effects on the flow field 
were evaluated with a viscous-inviscid interacting computational model. 
This model was constructed by coupling the explicit Euler solution 
procedure with a "lag-entrainment" boundary layer solution procedure in a 
global iteration scheme. The computational techniques have been used to 
compute the flow field for a long-duct turbofan engine nacelle at free- 
strea-n Mach numbers of 0.80 and 0.94 and angles of attack of 0° and 4°. 

The numerical experiments show that for predicting the flow inside 
the nacelle duct, the viscous effects are. extremely important and both the 
external and internal boundary layers and wakes must be simulated. The 
internal boundary layer creates a pressure gradient in the nacelle duct, 
the external boundary layer and wake displace the streamlines away from 
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the external surface of the nacelle, thereby reducing the compression at 
the trailing edge. This gives a less positive exit pressure, and hence a 
less .positive overall internal pressure level which agrees better with 
experiment than the inviscid computations. 

In contrast to the internal surface, viscous effects were relatively 
unimportant for predicting the flow on the external surface of the 
nacelle. Good agreement is shown between the computational results of 
both Euler numerical procedures and wind-tunnel .data on the external 
surface of the nacelle. The solutions at 0° angle of .attack exhibit the 
proper axisymmetric behavior even though the computational techniques are 
three dimensional. At 4° angle of attack., the solution has a. definite 
three-dimensional character. The calculations correctly predict the 
changing nature of the flow at the supercritical free stream Mach number 
of 0.94, and predict .the quantitative results at this Mach number very 
wel 1 . 

The solutions using both Euler computational techniques exhibited a 
total pressure loss on the internal surface of the nacelle. An 
investigation of the loss proved that it was the result o f the flow 
stagnating on the external surface and expanding around the sharp 
discontinuity in the surface of the nacelle at its leading edge. The 
studies indicate that the use of ... C-type grids could probably eliminate 
the loss, .although even with the C-grid, care.must be taken to use 
sufficient grid resolution to resolve the stagnation region near the 
leading edge. However, . for sharp leading edges or where the H-type grid 
Is otherwi se...desi rabl e , it appears that some (problem-dependent) total 
pressure loss is inherent to numerical Euler-equation solutions. 
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Finally, in comparing, the. two- numerical techniques, the explicit 

computational technique gives a more accurate solution for. the same number 

of grid points than the -implicit technique ._ The implicit computer code 

requires less computer storage while the explicit code runs faster. For 

the implicit computational procedure, it was found that a radiation. 

boundary condition imposed at the far-field and outflow boundaries gives 

better convergence of the computations than the condition p = p w . Even 

though reasonably accurate engineering solutions were obtained with the 

implicit computational procedure, a weak instability was discovered in it 

when applied to the. three-dimensional nacelle problem and the solutions 

never actually converged. This instability would severly restrict the use 

of this implicit method for studying problems of this type. The explicit 

computational technique is stable; however, its convergence rate is less i 

* 

than desirable. The weak instability of the implicit scheme coupled with . jj 

the slow rate of convergence of the explicit scheme provides motivation 

•i 

for further algorithm development. ! 

t ; 
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APPENDIX 


TWO-DIMENSIONAL INVESTIGATION OF THE' TOTAL PRESSURE LOSS 

In Chapter VI, part of the discussion centered upon a total pressure 
loss in the computed flow on the internal surface of the nacelle, a 
feature of the flow. that is physically incorrect. It was stated that the 
loss was first noticed for the solutions obtained with the alternating- 
di recti on -implicit computational procedure, and that at the time, it was 
attributed to incorrect treatment of the leading edge boundary. The 
explicit Runge-Kutta computational procedure does not contain a direct 
treatment of the leading edge boundary. Therefore, when solutions 
computed by the explicit procedure contained a similar loss in total 
pressure, an investigation to determine the cause of the loss was 
considered desirable. This appendix summarizes the resulting 
investigation. 

For the more detailed investigation into the source of the loss, an 
unswept, untapered wing, was chosen as the test configuration. Using the 
unswept, untapered wing resulted in a physical problem with two- 
dimensional flow, but retained the three-dimensionality of the computa- 
tional procedure. Symmetry of the boundary conditions could be 
maintained. Also, by adopting a symmetrical airfoil, and by comparing the 
solutions on the top and bottom surfaces of the wing, the. coding of the 
surface boundary condition inside the nacelle Could be verified. This was 
possible because the boundary conditions, and hence the computer codings, 
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for the bottom surface of the wing and for the Inside surface of tne 
hace'i le are Identical . 

Solutions for the wing were obtained using the explicit Runge-Kutta 
algorithm for both a. symmetrical airfoil and a cambered airfoi-1 with H- 
grids and C-grids. Figure A-l shows examples of the two types of grids. 
The H-grid resembled as closely as possible the grid used for the 
nacelle. All calculations presented in the appendix are at an angle. of 
attack of 0.0°. 

Figure A-2 presents the solution, at a free, stream Mach number of 
0.80 for a wing with a symmetrical airfoil and an H-grid. The. solution 
is symmetrical about the cord line, and, even though there are regions of 
supersonic flow, the total pressures on both surfaces are essentially at 
the free stream value. The slight variation in the total pressures where 
the flow compresses from the negative pressure peak, is possibly the result 
of excessive numerical dissipation.. The symmetry of this solution, and 
the equality of the total pressures on the surface to the free stream 
total pressure strongly indicates that the surface boundary condition is 
programed correctly. 

The computed total pressures for the wing with a cambered airfoil 
having the same cross-section as the nacelle and an H-grid are presented 
in figure A-3. This solution is also at a free stream Mach number of 
0.80 and an angle of attack of 0.0°. The figure shows that there is a 
larger total pressure loss on the bottom surface of the wing than on the 
corresponding inside surface of the nacelle, and there is even a. -loss on 
the top surface. 

In order to investigate the effect of rounding the leading edge, 
it's contour was changed to a circular arc which was tangent to the top 








2 



Figure A-3. Computed total pressures for the wing with a 
cambered airfoil and an H-grid. (Kl = 0.80, 
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and bottom surfaces of the airfoil, and. a C-grid was adopted. Figure A- 
4(a) presents the total pressure distributions for this configuration at 
the same free stream conditions as the previous calculations. It shows 
that while the the total pressure is considerably closer to the free, 
stream value than for the H-grid, it is now too high. Whether the calcu- 
lated total pressures on the surface are higher or lower than the free 
stream value depends upon the interaction between the implementation of 
the boundary conditions and the numerical dissipatiorr'- . 

An examination of the Mach number contours for this calculation, 
presented in figure A-4(b) and (c), reveals two interesting features of 
the flow. The first feature is that the cambered airfoil is at 0.0° angle 
of attack, and yet the flaw stagnates on the upper surface of the wing 
near the leading edge. In addition, as the flow expands around the 
leading, edge, a supersonic bubble with relatively high Mach numbers is 
created on the bottom surface. The gradients in the supersonic bubble 
coupled with the grid resolution and the numerical dissipation probably 
produce the change in total pressure on the bottom surface of the wing. 
Possibly a finer grid in this region could alleviate the problem. 

However, the present calculations resulted in near saturation of the 
incore storage of the Cyber 203 computer. 

By eliminating the severe expansion around the leading edge, the 
total pressure losses on the surface of the wing can be prevented. This 
has been demonstrated for a symmetrical airfoil using the H-grid and was 
presented In figure A-2. It will also be demonstrated, for the C-grid in 
two ways, first by going to a symmetrical airfoil as in the case for the 
H-grid, and Second, by lowering the free stream Mach number. Figure A-5 
presents the solution for the wing with the symmetrical airfoil and the C- 
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L CONTOUR =0.10 

M = 1.00 



(c) Detail of the Mach number contours at the leading edge 

Figure A-4. Concluded. 
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(c) Total pressure distributions. 


Figure A-5. Concluded 


grid at a free stream Mach number of 0.80.. The pressure coefficients, 
which are presented in. part (a) of the figure, agree well with the corre- 
sponding pressure distributions calculated with the H-grid. Figure A- 
5(b), which shows the Mach number contours, illustrates that, the 
stagnation point is at the forward most point of the airfoil, and that the 
strength of the expansion at the leading edge has been reduced. ... The 
resulting total pressure distributions are presented in part (c) of the 
figure, and. show that the total pressure on the surface of the wing 
remains close to the free stream value. 

By reducing the free-stream Mach number, it is also possible to 
reduce the streogtjh of the expansion. The results for a free stream Mach 
number of 0.40 are presented in figure A-6. The pressure coefficients, 
presented in pa.rt (a), look reasonable although there are no experimental 
data with which to compare them. .The Mach number contours in the vicinity 
of the leading edge, and the total pressure distributions, parts (b) and 
(c) respectively, illustrate the reduced strength of the expansion at the 
leading edge, and the resulting surface total pressures which are 
essentially at the free stream value. 

An investigation to determine the cause of the total pressure loss 
in the computed flow on the internal surface of the nacelle has been 
summarized in . this appendix. It was found that the loss was associated 
with the flow expanding around the sharp leading edge of the nacelle and 
thus creating severe local gradients in the flow field. It was determined 
that by reducing the severity of these gradients, the total pressure loss 
could be prevented. Although the specific approaches used for the wing 
(C-grid and reduced Mach number) were not attempted for the nacelle, an 
alternate approach in which the nacelle geometry was modified such that 









stagnation occured precisely at the leading edge was tried* This jnodl fled 
nacelle produced similar results in that the loss In total pressure was 
eliminated (see.Cnapter VI). The implication is that almost any modifica- 
tion, to the leading -edge treatment that either reduces or resolves the 
flow gradients should produce similar results. 
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